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1  Introduction 


In  recent  publications  [2,  3,  4,  5] a  new  globally  convergent  numerical  method  for  a  Coefficient 
Inverse  Problems  (CIP)  for  a  hyperbolic  Partial  Differential  Equation  (PDE)  was  developed 
analytically  and  tested  numerically  on  computationally  simulated  data.  Thus,  the  goal  of 
the  effort  described  this  grant  was  to  verify  the  performance  of  that  technique  on  a  set  of 
experimental  data  measured  in  the  picosecond  scale  regime.  The  Co. -PI  and  his  graduate 
student  have  obtained  measurements  of  tomographic  time  resolved  data.  The  total  timing 
of  measurements  per  one  detector  location  was  12.3  nano-seconds  per  one  detector  location. 
Recall  that  1  pico-second  (ps)=10-12  second=lCT3  nano-seconds  (ns).  The  measured  data 
were  pre-processed  using  a  radically  new  data  processing  procedure.  The  processed  data 
were  used  then  as  Dirichlet  boundary  conditions  for  elliptic  PDEs  derived  in  [2,  3].  Next, 
images  were  obtained  by  the  algorithm  of  [2,  3]. 

Remark.  Because  of  some  technical  difficulties,  we  have  “double  numbering”  for  some 
citations  below:  numbers  [1]- [4]  in  subsection  12.1  are  the  same  as  [l]-[4]  in  subsection  12.2, 
although  they  really  denote  different  citations.  Nevertheless,  it  is  always  clear  from  the 
content  “what  is  what” . 

1.1  Our  results  in  brief 

Only  semi-blind  data  were  used.  Namely,  it  happen  that  we  knew  locations  of  dielectric 
inclusions,  although  this  information  is  not  used  in  our  algorithm.  However,  we  did  not 
know  values  of  refractive  indexes  of  those  inclusions.  Therefore,  we  use  the  term  blind, 
everywhere  below  when  talking  about  refractive  indexes  of  inclusions  used  in  experiments. 
After  computational  results  were  obtained,  refractive  indexes  of  inclusions  were  measured  a 
posteriori  by  two  well  established  methods.  Comparison  of  measured  and  blindly  computed 
refractive  indexes  has  consistently  demonstrated  a  very  good  accuracy  of  computed  ones.  It 
is  an  opinion  of  the  authors  that  the  latter  fully  verifies  the  technique  of  [2,  3,  4,  5]. 

We  have  obtained  quantifiable  non-destructive  images  of  dielectric  inclusions  hidden  in 
otherwise  slowly  changing  background.  “Quantifiable”  means  that  the  values  of  refractive 
indexes  in  those  dielectric  inclusions  are  accurately  imaged.  The  accurate  imaging  of  these 
values  is  an  important  ingredient  in  the  goal  of  the  identification  of  those  abnormalities. 
Indeed,  these  values  might  help  to  differentiate  between  various  types  of  dielectrics.  Po¬ 
tential  applications  of  the  quantifiable  imaging  of  dielectrics  are  in  checking  out  baggages 
in  airports,  a  stand-off  detection  of  potential  explosives  hidden  under  clothing,  imaging 
of  antipersonnel  land  mines,  etc..  Indeed,  it  is  well  known  from,  e.g.  tables  presented 
at  http://www.clippercontrols.com/info/dielectric_constants.html^  1  that  explosives  have 
much  higher  dielectric  constants  than  ones  of  regular  materials.  The  framework  of  the  tech¬ 
nique  of  [2,  3]  can  be  extended  to  the  case  when  both  the  dielectric  permittivity  and  the 
electric  conductivity  coefficients  are  unknown,  although  we  have  not  worked  out  specific 
details  yet. 

The  publication  with  experimental  results  was  highly  rated  by  referees  and  became  a 
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featured  article  of  Inverse  Problems  [2].  In  addition  to  the  main  result  of  [2],  we  also  describe 
here  three  directions  of  our  ongoing  effort  supported  by  this  funding.  The  ultimate  goal  of 
this  effort  is  to  refine  results  of  [2]  in  terms  of  refining  shapes  of  imaged  inclusions. 

On  the  first  stage  of  this  project,  the  above  blind  images  have  shown  locations  of  inclusions 
and  values  of  refractive  indexes  in  them.  In  this  subsection  we  describe  our  effort  to  get 
better  sizes  of  inclusions  from  experimental  data.  To  do  this,  we  have  applied  the  two-stage 
numerical  procedure,  which  has  been  developed  under  funding  of  W911NF-08- 1-0470  during 
the  year  2009  [?].  So,  now  this  development  pays  off.  On  the  first  stage  of  this  procedure 
the  globally  convergent  method  of  [2]  provides  a  guaranteed  good  first  approximation  for  the 
solution.  And  on  the  second  stage  this  approximation  is  refined  via  a  subsequent  application 
of  a  locally  convergent  Adaptive  Finite  Element  technique. 

1.2  Comparison  with  conventional  approaches 

In  2005  and  2009  Inverse  Problems  has  published  two  special  issues  devoted  to  reconstruc¬ 
tions  of  both  dielectric  and  conductive  abnormalities  from  experimental  data  [7,  13,  22,  23]. 
While  the  data  of  [22]  of  2005  were  designed  for  two-dimensional  images,  the  data  [23]  of  2009 
were  collected  for  three-dimensional  ones.  These  experimental  data  were  provided  by  Fresnel 
Institute  (Marseille,  France).  A  number  of  good  quality  imaging  results  was  published  in 
these  issues.  We  now  list  three  main  differences  between  our  work  and  ones  in  [7,  13,  22,  23]. 
These  differences  emphasize  the  advantage  of  our  approach  over  existing  ones. 

1.  Sources/Detectors  Configurations.  In  our  case  only  a  single  location  of  the 
source  is  used,  and  tomographic  measurements  are  performed  only  on  one  side  of  a  prism, 
which  is  opposite  to  the  source  location,  see  Figure  1  in  section  5.  Our  source/detectors 
configuration  uses  the  minimal  amount  of  information.  In  addition,  our  source/detectors 
configuration  can  be  easily  transformed  in  the  most  interesting  case  (from  the  practical 
standpoint)  of  backscattering  data,  i.e.  to  the  case  of  stand-off  detection.  The  technique  of 
[2,  3,  4,  5]  can  also  be  extended  to  this  case,  see  subsection  9.2.  In  [22,  23]  the  source  is 
moved  all  around  the  medium  of  interest  and  the  data  are  collected  all  around  a  circle.  The 
latter  is  of  course  inconvenient  for  Army  applications. 

2.  Convergence  Analysis.  The  global  convergence  of  our  algorithm  is  rigorously  guar¬ 
anteed  [2,  3].  We  work  with  a  fully  nonlinear  problem  and  do  not  use  neither  a  linearization, 
nor  an  assumption  that  the  starting  point  is  located  close  to  the  correct  solution,  nor  an 
assumption  about  a  priori  knowledge  of  the  background  medium.  The  only  a  priori  knowl¬ 
edge  we  use  is  that  the  medium  outside  of  our  domain  of  interest  is  air  and  that  the  relative 
dielectric  permittivity  in  the  domain  of  our  interest  is  not  less  than  the  one  in  the  air,  i.e. 
that  the  EM  wave  propagates  in  the  domain  of  interest  slower  than  in  the  air.  Contrary  to 
this,  works  of  [22,  23]  use  some  variations  of  the  small  perturbation  approach,  e.g.,  Newton¬ 
like  techniques.  Convergence  of  the  small  perturbation  approach  cane  be  proven  only  under 
the  condition  that  the  starting  point  is  located  close  to  the  correct  solution. 

3.  Data  Collection.  We  have  measured  the  time  resolved  data  generated  by  a  100  ps 
electric  pulse.  In  [7,  13,  22,  23]  the  data  are  collected  on  several  frequencies  by  the  so-called 
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Vector  Network  Analyzer  technique.  While  formally  measurements  on  several  frequencies 
are  equivalent  with  time  resolved  ones  via  the  Fourier  transform,  it  is  the  time  resolved  signal 
which  provides  the  ultimate  result. 

1.3  Arrangements 

Since  this  is  only  a  nine  months  funding,  the  Co-PIs  have  realized  prior  to  the  starting 
date  of  funding  that  it  would  be  highly  risky  to  count  that  a  new  research  software  for  the 
globally  convergent  method  would  be  developed  and  thoroughly  tested  during  this  period. 
In  addition,  because  of  funding  limitations  ($50,000),  we  have  realized  that  it  is  impossible 
to  buy  a  very  expensive  equipment  we  needed.  Given  these  timing  and  financial  constraints, 
it  is  unlikely  that  the  project  would  succeed  without  two  key  arrangements  described  in  this 
subsection. 

The  first  arrangement  was  with  the  research  software  for  the  method  of  [2,  3].  Fortunately, 
Dr.  L.  Beilina  (Chalmers  University,  Gothenburg,  Sweden),  who  has  been  collaborating  with 
the  PI  since  2007  (with  joint  publications  cited  as  [2,3]  in  subsection  12.1  and  as  [2-5]  in 
subsection  12.2)  has  kindly  agreed  to  let  us  to  use  her  research  software,  which  she  has  been 
developing  for  a  number  of  years.  So,  the  Postdoctoral  Research  Associate  of  the  PI,  Dr. 
Natee  Pantong  has  visited  Beilina  for  one  month  of  September  2009,  using  travel  funds  from 
W911NF-08-1-0470.  The  goal  of  his  visit  was  to  learn  how  to  work  with  this  software.  So, 
results  described  below  in  sections  6  and  7  were  obtained  using  that  software  as  the  main  part 
of  the  computational  work.  In  addition,  Dr.  Pantong  has  written  a  supplemental  software 
for  experimental  data  processing. 

The  second  crucial  arrangement  was  with  the  equipment.  We  have  decided  to  rent  it  for 
one-two  months  period  rather  than  buy  it. 

2  Forward  and  Inverse  Problems 

As  the  forward  problem,  we  consider  the  following  Cauchy  problem 

£r(x)utt  =  A u,  in  M3  x  (0,oo) ,  (1) 

u  (x,  0)  =  0,  ut  (x,  0)  =  5  (x  —  x0)  •  (2) 

Here  er(x)  is  the  spatially  variable  dielectric  constant  (relative  dielectric  permittivity), 

£r(x)  =  y/er{x )  =n(x)  =  (3) 

£0  C  (X) 

where  £0  is  the  dielectric  permittivity  of  the  vacuum  (which  we  assume  to  be  the  same  as  one 
in  the  air),  £  (x)  is  the  spatially  variable  dielectric  permittivity  of  the  medium  of  interest, 
n  (x)  is  the  spatially  variable  refractive  index  of  the  medium  of  interest,  c  (x)  is  the  speed  of 
the  propagation  of  the  EM  field  in  this  medium,  and  Cq  is  the  speed  of  light  in  the  vacuum, 
which  we  assume  to  be  the  same  as  one  in  the  air.  Thus,  the  refractive  index  shows  how 
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slower  the  EM  field  propagates  in  the  medium  of  interest  compared  with  the  air.  We  assume 
below  that  n  (x)  >  1. 

We  point  out  that  equation  (1)  can  be  derived  from  the  Maxwell’s  system  only  in  two 
cases:  (a)  if  er  =  const.  >  0  and  (b)  in  the  2-D  case  [8].  However,  neither  of  these  two  is  in 
place  in  our  experiments.  Thus,  we  call  (1)  a  simplified  mathematical  model  of  our  process. 

Let  11  C  M3  be  a  convex  bounded  domain  with  the  boundary  dfl  G  C3.  We  assume  that 
the  coefficient  er  ( x )  of  equation  (1)  is  such  that 

er  ( x )  >  1,  er  (x)  =  1  for  x  G  M3\fi,  (4) 

£r  ( x )  G  C2  (M3)  .  (5) 

It  is  well  known  that  the  question  of  the  existence  of  the  fundamental  solution  (1),  (2)  is 
a  very  challenging  one.  Regardless  on  an  extensive  effort  in  the  past  of  such  people  as  J. 
Hadamard,  S.L.  Sobolev  and  V.G.  Romanov,  this  question  is  currently  positively  addressed 
only  under  the  condition  that  the  coefficient  £r  ( x )  is  sufficiently  smooth  and  geodesic  lines 
generated  by  this  coefficient  are  regular,  see  works  of  Romanov  [20,  21]  and  citations  of 
books  J.  Hadamard  and  S.L.  Sobolev  in  them.  Hence,  it  is  inevitable  that  in  the  theoretical 
derivations  of  [2,  3]  both  (5)  and  the  regularity  of  geodesic  lines  were  assumed.  Note  that 
we  need  these  assumptions  only  to  make  sure  that  the  solution  of  the  problem  (1),  (2)  exists 
and  that  its  Laplace  transform  (7)  has  a  certain  asymptotic  behavior,  see  (16).  However, 
we  actually  do  not  use  these  assumptions  in  our  computational  practice  and  verify  (16) 
computationally,  see  subsection  7.2  in  [2],  Likewise,  although  in  our  experiments  the  function 
£r  ( x )  has  a  discontinuity  at  the  boundary  of  a  dielectric  inclusion,  which  contradicts  to  (5), 
our  reconstruction  method  still  works. 

Inverse  Problem.  Suppose  that  the  coefficient  £r(x)  satisfies  (4)  and  (5)  and  cor¬ 
responding  geodesic  lines  are  regular.  Assume  that  the  function  £r  ( x )  is  unknown  in  the 
domain  fL  Determine  the  function  £r  ( x )  for  x  G  fl,  assuming  that  the  following  function 
g  ( x ,  t )  is  known  for  a  single  source  position  x0  ^ 

u  (x,  t)  =  g  ( x ,  t) ,  V  (x,  t )  G  dfl  x  (0,  oo) .  (6) 

In  our  application  the  assumption  £r  ( x )  =  1  for  x  G  M3\h2  means  that  one  has  air 
outside  of  the  medium  of  interest  fh  The  inequality  £r  ( x )  >  1  is  because  is  that  one  should 
bound  the  coefficient  £r  ( x )  from  the  below  by  a  positive  number  to  ensure  that  the  operator 
in  (1)  is  a  hyperbolic  one  on  all  iterations  of  our  numerical  procedure.  In  addition,  this 
inequality  means  that  the  EM  wave  propagates  slower  in  the  domain  of  interest,  compared 
with  the  air,  see  (3).  The  function  g  (x,  t )  models  time  dependent  measurements  of  the  wave 
field  at  the  boundary  of  the  domain  of  interest. 

3  Outline  of  the  Globally  Convergent  Numerical  Method 

We  outline  in  this  section  the  method  of  [2,  3]  as  well  as  the  global  convergence  theorem  of 
[3],  which  is  more  advanced  than  one  in  [2],  We  refer  to  these  publications  for  details. 
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3.1  Outline 

Consider  the  Laplace  transform  of  the  functions  u, 


OO 


e  stdt  :=  jC(u)  ,  for  s  >  s  =  const.  >  0, 


(7) 


o 


where  s  is  a  sufficiently  large  number  such  that  the  integral  (7)  converges  together  with 
corresponding  (x,  f)-derivatives.  We  call  the  parameter  s  pseudo  frequency.  Note  that  we  do 
not  use  the  inverse  Laplace  transform  in  our  method,  since  approximations  for  the  unknown 
coefficient  are  obtained  in  the  pseudo  frequency  domain.  We  obtain  from  (1),  (2) 


Aw  —  s2c  (x)  w  =  —5(x  —  x0) , 
lim  w  (x,  s )  =  0. 

|  x  |  — >oo 


(8) 

(9) 


To  prove  (9),  one  can  first  consider  a  Laplace-like  transform  of  the  function  u(x,t)  [?], 


o 

It  was  shown  in  [?]  that  the  function  u  satisfies 

er  (x)  ut  =  A u,  u  (x,  0)  =  S  (x  —  Xo)  • 

Next,  for  sufficiently  large  s 


(10) 


OO 


(11) 


o 


It  follows  from  the  classic  estimate  (6.13)  of  Chapter  4  of  [?]  for  the  fundamental  solution  of 


the  parabolic  PDE  that  the  function  \u\  can  be  estimated  from  the  above  via  the  solution 


of  another  parabolic  equation  with  constant  coefficients.  Finally,  the  Laplace  transform  (11) 
for  the  latter  solution  can  be  calculated  via  an  explicit  formula,  and  this  formula  implies 
decay  as  |x|  —>  oo. 


It  follows  from  the  classic  theory  of  PDEs  that  for  every  s  >  s  there  exists  unique  solution 


w  G  C2+a  (|x  —  Xo|  >  7) ,  V7  >  0.  Here  and  below  Ck+1,  a  G  (0, 1) ,  k  >  0,  an  integer,  are 
Holder  spaces  [18].  Since  the  fundamental  solution  u  of  the  problem  (10)  is  positive  (Theorem 
11  of  Chapter  2  of  [11]),  then  by  (11)  w(x,s )  >  0  for  sufficiently  large  s  in  (7).  Hence,  we 
can  consider  the  function  v(x,s)  =  s~2  lnu>(x,  s),  which  is  the  Liouville  transform  of  the 
function  w.  The  function  v  satisfies  the  following  conditions 


Av  +  s2  ( Vu )2  =  £r  (x)  ,  X  G  n, 
V  |  dn  =  i>{x,s), 


(12) 
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where  the  function  i/j  is  generated  by  the  function  g  in  (6).  While  (12)  was  obtained  using  the 
Liouville  transform,  a  new  step  was  proposed  in  [2].  The  idea  of  this  step  has  a  root  in  the 
idea  of  applications  of  Carleman  estimates  to  proofs  of  uniqueness  theorems  for  CIPs  [15,  16]. 
Namely,  we  eliminate  the  unknown  coefficient  er  (x)  from  equation  (12)  via  differentiating 
this  equation  with  respect  to  s.  Let  q  (x,  s )  =  dsv  (x,  s )  .  Then  we  obtain  from  (12)  that  the 
function  q  (x,  s )  satisfies  the  following  nonlinear  integral  differential  PDE  containing  Volterra 
integrals  with  respect  to  s, 

S 

Ag  —  2s2Vg  ■  j  Vg  (x,  r)  dr  +  2s 

S 

}  (13) 

+  2 s2  Vg  •  W  —  2 s W  •  /  Vg  (x,  r)  dr  +  2s  (VP)2  =  0,i  6  O, 

S 

q  |n=  Op  s)  ,  (x,  s)  E  dfl  x  [s,  s] , 


S 

j  Vg  (x,  r)  dr 

S 


where  ^  (x,  s)  =  dsij)  (x,  s ) .  Here  s  is  the  truncation  pseudo  frequency  of  integrals,  which 
serves  as  one  of  regularization  parameters  of  our  method.  Still,  instead  of  just  truncating  the 
integral  via  setting  its  complement  to  zero,  we  use  the  function  V  (x,  s) ,  which  complements 
the  rest  of  the  integral,  i.e. 


S 


=  -Jq(*,T)dr  +  V(X,s), 

(14) 

s 

T  r  /  — \  /  _x  In  w(x,s) 

V  [X,  s)  =  V(X,s)  —  - ~2 - . 

(15) 

We  call  V  (x,  s)  the  “tail  function” ,  and  it  is  unknown.  Hence,  equation  (13)  has  two  unknown 
functions,  q  and  V.  The  reason  why  we  can  approximate  well  both  of  them  is  that  we  treat 
them  differently.  While  we  approximate  the  function  q  from  inner  iterations,  the  function  V 
is  approximated  via  outer  iterations.  In  fact,  numerical  solution  of  the  problem  (13)  is  the 
most  challenging  issue  in  this  method.  It  follows  from  (15)  and  Lemma  2.1  of  [2],  which  is 
actually  based  on  Theorem  4.1  of  [20]  as  well  as  on  the  work  [21]  of  the  same  author,  that 


_d_ 

dsk 


V  (x,  s ) 


C2+“(o) 


0,1,... 


(16) 


The  problem  (13)  is  solved  via  a  layer  stripping  procedure  with  respect  to  the  pseudo 
frequency  s.  Consider  partition  of  the  interval  [s,  s]  into  N  small  subintervals  of  the  width 
h  =  sn_!  —  sn,  where  s  =  sjy  <  sjv-i  <  •••  <  Si  =  s.  Assume  that  q(x,s)  is  a  piecewise 
constant  function  with  respect  to  s,  q  (x,  s )  =  qn  (x)  for  x  E  (sn,  sn_i] .  Consider  the  Carleman 
Weight  Function  (CWF)  e^s~Sn~i\  where  /i  >>  1  is  a  large  parameter  which  should  be 
chosen  in  computations.  We  multiply  equation  (13)  by  this  function  and  integrate  with 
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respect  to  s  E  (sn,  sn_ i) .  Hence,  we  obtain  the  Dirichlet  boundary  value  problem  for  the 
following  coupled  system  of  nonlinear  elliptic  PDEs  of  the  second  order 


71—1 


Ln  (On)  •  ^ Qn  l,n  v  ^ Qi  j  ^ Qn  ^  % Qn  Bn  Qn) 

- A2,nh2  ( Vft  (*))  +  2A2?n Vy  (h  ]r  W)  -  A2,n  (' W)2  , 


,  2=1 


^n- 1 


(17) 


2=1 


an  =  ^n(,x):=-  /  i/;(x,s)ds,  n  = 


where  the  function  ^{x,s)  is  taken  from  (13).  In  (17)  A^n,  A2,n,  Bn  are  certain  numbers 
depending  on  /i,  h,  n  and  x  >  0  is  a  small  parameter  of  ones  choice.  This  parameter  is 
introduced  to  obtain  a  better  stability  of  the  problem  (17)  because  of  the  maximum  principle, 
see  §1  in  Chapter  3  of  [18].  It  is  important  that  lim^,^  Bn  =  0  uniformly  for  all  n  due  to  the 
presence  of  the  CWF.  Hence,  the  presence  of  the  CWF  with  /i  >>  1  mitigates  the  influence 
of  the  nonlinear  term  (Vgn)2  ,  which  enables  us  to  solve  the  boundary  value  problem  for 
each  qn  iteratively  via  solving  a  linear  elliptic  problem  on  each  step.  Still,  the  computational 
experience  shows  that  we  cannot  take  /i  exceedingly  large,  which  would  effectively  turn 
equations  (17)  into  linear  ones.  Suppose  that  we  have  approximated  the  function  qn  (x) . 
Then  we  find  the  approximation  c”  (re)  for  the  function  er  (#)  via  backwards  calculation 
using  (12)  as 

Jn)  (  X  _  /  fn  (x)  :=  Avn  +  s2n  (Vvnf  ,  x  E  Q,  if  fn  (x)  >  1, 
r  {  )  \  1  if  fn  (x)  <  1, 

n 

Vn  (x)  =  -h  ^2  Qi  (X)  +  Vn  (x)  , 

2=1 

where  Vn  (x)  is  the  corresponding  approximation  for  the  tail  function.  We  make  the  cut-off 
to  unity  in  (18a)  because  of  (4). 

Since  the  equation  for  each  function  qn  in  (17)  depends  only  on  functions  q±,  ...,qn,  then 
these  elliptic  Dirichlet  boundary  value  problems  can  be  solved  sequentially:  first  one  should 
approximate  q\ ,  next  approximate  q2,  etc..  We  have  inner  and  outer  iterations  to  solve  these 
problems.  While  functions  qn  are  found  from  inner  iterations,  approximations  for  the  tail 
function  V  are  found  from  outer  iterations.  First,  we  choose  the  first  approximation  VW  (x,  s) 
for  the  tail  function.  In  [2]  Vlti  =  0  was  chosen.  While  we  can  still  do  so,  we  have  discovered 
that  the  process  converges  faster  if  we  choose  the  approximation,  which  corresponds  to  the 
case  of  the  known  value  of  the  function  er  =  1  outside  of  the  domain  fh  So,  let  w  (x,s)  be 
the  solution  of  the  problem  (8),  (9)  with  er  =  1.  Using  (15),  we  set  [3] 


(18a) 

(18b) 


To  approximate  the  function  gi,  we  first  iterate  with  respect  to  the  nonlinear  term  (Vgi)2 
in  (17)  and  find  functions  q\  (x)  this  way  setting  in  (17)  V  :=  V^i,  where  Vj  j  is  defined  in 
(19).  So,  we  solve  the  Dirichlet  boundary  value  problem  (17)  for  functions  q\  (x)  via  setting 
for  the  nonlinear  term  (Vgn)2  :=  (Vg^1)  ,  q®  —  0.  We  iterate  with  respect  to  the  nonlinear 
term  until  convergence  occurs.  The  resulting  function  is  denoted  as  gi.i. 

For  n  =  1,  we  do  not  iterate  with  respect  to  the  nonlinear  term  anymore.  Instead,  we 
iterate  with  respect  to  the  tail  as  follows.  Suppose  we  have  obtained  the  pair  (qlti,  Vj^)  • 
Then  we  find  the  approximation  si1’^  ( x )  for  the  target  coefficient  er  ( x )  via  the  backwards 
calculation  (18a, b)  with  the  obvious  replacement  of  indexes.  Next,  we  solve  the  problem  (1), 
(2)  with  er  :=  £r  't'> ,  calculate  the  Laplace  transform  W\^+i  (x,  s)  of  its  solution  and  set  the 
new  approximation  for  the  tail  as  V\,l+\  ( x )  :=  s~2Wi,i+i  ( x,s ) .  Next,  we  solve  the  boundary 
value  problem  (17)  with  V  :=  Vlji+ 1,  (Vgra)2  :=  (Vgp;)2  and  obtain  the  function  g1)i+1  this 
way.  We  continue  this  process  until  convergence  occurs.  Suppose  that  convergence  occurs  at 
i  :=  m\.  Then  we  set  ^gi,  £^\  :=  ^gi>mi,  ,  for  the  nonlinear  term  in  (17) 

we  set  (Vgn)2  :=  (Vg2)2  :=  (Vgi)2  .  To  find  the  function  g2,  we  repeat  the  above  process  for 
n  =  2,  etc.,  until  convergence  occurs.  So,  for  each  n  we  iterate  with  respect  to  the  nonlinear 
term  only  to  approximate  gnjl  and  then  we  iterate  with  respect  to  tails.  We  can  rigorously 
prove  convergence  of  iterations  with  respect  to  the  nonlinear  term  for  each  n  (Theorem  1). 
However,  we  cannot  prove  convergence  with  respect  to  tails.  Still,  we  have  observed  the 
latter  numerically  and  our  numerical  convergence  criteria  are  specified  in  subsection  7.1. 

Since  we  have  defined  functions  £r'J)  in  (18a, b)  only  for  x  G  0,  we  explain  now  how  we 
can  extend  these  functions  in  M3\fl  in  such  a  way  that  the  resulting  function  belongs  to 
Ca  (M3) .  We  need  this  explanation  only  for  the  sake  of  our  global  convergence  theorem  in 
the  next  subsection,  since  in  the  computational  practice  we  just  extend  these  functions  as 
(x)  :=  1  for  x  G  M3\f2.  Let  fT  C  be  a  convex  subdomain  such  that  dist  (dfl1,  dfl)  >  0, 
where  dist  (dfl1,  dfl)  is  the  distance  between  these  boundaries  in  the  Hausdorf  sense.  We 
assume  that  the  number  dist  (<9fT ,  dfl)  is  rather  small,  i.e.  «  fl.  It  is  well  known  from 
the  Real  Analysis  course  that  one  can  choose  such  a  function  x  ix)  £  C°°  (M3)  that 

(  1  in  O', 

X  (x)  =  <  between  0  and  1  in  Q\Q', 

I  0  in  M3\f2. 


So,  suppose  that  in  (18a, b)  the  function  (x)  G  Ca  (fl)  .  We  define  the  function 


as  Er  [X) 

1  for  x  €  I 


=  (1  —  X  (^))+X  (x)  (a^) ,  Vx  G  M3.  Then  ein'^ 

x)  in  f2'. 


G  C° 


),£(rn’i)>l ,3""  (X)  = 


^3\f2  and 


\X)  =  er”’0 


x) 

x)  = 


3.2  The  global  convergence  theorem 

To  formulate  this  theorem,  we  need  to  introduce  the  definition  of  the  exact  solution  first.  We 
assume  that  there  exists  a  coefficient  e*  ( x )  satisfying  condition  (4),  (5),  and  this  function 
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is  an  exact  solution  of  our  Inverse  Problem  with  the  “ideal”  errorless  data  g*(x,t)  in  (6). 
The  Laplace  transform  (7)  of  the  function  g*  (x,t)  leads  to  the  exact  function  p*  (x,s)  = 
w*  (x,  s) ,  V  (x,  s)  G  dtt  x  [s,  s] .  Let 


\  9 
q  m  =  Ts 


In  [w*  (x,  s)] 


Hence,  V*  (x,  s )  is  the  exact  tail  function.  The  function  q*  satisfies  an  obvious  analogue  of 
equation  (13)  with  the  following  boundary  condition 


q*  (x,  s )  =  ijj*  (x,  s ) 


1  dp* 
p*s2  ds 


2  In  p* 


,  (x,  s )  G  dQ  x  [s,  s] . 


It  easily  follows  from  the  above  that  the  function  q*  (x,  s )  G  C2+a  (H)  x  C°°  [s,  s]  .  First,  we 
approximate  functions  q*  (x,  s)  and  ifj *  (x,  s )  via  piecewise  constant  functions  with  respect  to 
s  G  [s,  s]  as 

Sn  —  1  $n  —  1 

qn(x)  =  \  J  Q*  (X’  S)  dsi  (X)  =  \  J  ^ *(x,s)ds .  (21) 

Sn  Sn 

Hence,  for  n  =  1,  ...,iV;s  G  [sn,sn_i]  we  have  q*(x,s)  =  g*  (x)  +  Qn  (x,  s) ,  -0*  (x,  s)  = 
V’n  (a;)+^/fi  (x> s ) )  where  functions  Qn,  are  such  that  | Qn  (x,  s)|2+Q  <  C*h,  (x,  s)|2+Q,  < 

C*h,  for  s  G  [sn,  sn_i] ,  where  the  constant  C*  =  C*  ^||g*||C2+a(n)xci[s  s])  >  0  depends  only 
on  the  C2+a  (hi)  x  Cl  [s,  s]  norm  of  the  function  q*  (x,s).  Here  and  below  |-|fc+Q  denotes  the 
norm  in  the  space  Ck+a  (hi)  .  We  can  assume  that 


max 

l<n<iV 


2+«  <  c* 


(22) 


and  without  a  loss  of  generality,  we  assume  that 


C*  >  1. 


(23) 


By  the  Tikhonov  concept  [24],  the  constant  C*  should  be  known  a  priori.  It  is  reasonable 
to  assume  that  C*  is  independent  on  s,  although  we  do  not  use  this  assumption.  By  (21) 
g*  ( x )  =  4>n  (x)  d  £  dVL.  We  assume  that  the  function  g(x,t )  in  (6)  is  given  with  an  error. 
This  naturally  produces  an  error  in  the  function  (x,  s)  in  (13).  An  additional  error  is 
introduced  due  to  the  averaging  in  (21).  Hence,  it  is  reasonable  to  assume  that 


^ n  0)  -  ^ n  0) 


<C*{a  +  h), 

C2+“(90) 


(24) 


where  a  >  0  is  a  small  parameter  characterizing  the  level  of  the  error  in  the  data  (x,  s ) . 
We  assume  that 


s  >  1,  nh>  1, 


(25) 
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(26) 


and  introduce  the  positive  constant  M*  =  M*  ^ ||g*|| c2+a ^xc^ss]  >  SJ  =  (C*,s)  by 

M*  =  16  C*s2. 


Consider  the  Dirichlet  boundary  value  problem  in  the  domain  hi 

3 

A u  +  ^2  bj(x)uXj  -  d(x)u  =  pi  ( x ) ,  u  |9n=  P2  (x)  E  C2+a  (dFl) . 
3= 1 


Assume  that  functions  bj,d,pi  E  Ca  (hi)  ,d(x)  >  0;  max  (\bj\a,  |d|Q)  <  1.  By  the  Schauder 
theorem  [18],  there  exists  unique  solution  u  E  C2+a  (hi)  of  this  boundary  value  problem,  and 
with  a  constant  K  =  K  (fl)  >  1,  depending  only  on  the  domain  fl,  the  following  estimate 
holds 


u 


2+a 


<  K 


\'Pl\a  +  \\p2  II C,2+“(9Q) 


(27) 


Theorem  1  (global  convergence)  [3].  Let  Q  C  R3  be  a  convex  bounded  domain  with  the 
boundary  dLl  E  C3.  Suppose  that  inequalities  (22)- (25)  hold.  Let  the  exact  coefficient  e*  (x) 
satisfies  conditions  (4),  (5).  For  any  function  er(x)  E  Ca  (M3)  such  that  er(x)  >  1  in  Ll 
and  £r(x)  =  1  in  M3\h2  consider  the  solution  u£r(x,t )  of  the  Cauchy  problem  (1),  (2). 
Let  w£r  (x,s)  E  C2+a  ({|a:  —  xfi  >  7}) ,  V7  >  0  be  the  Laplace  transform  (7)  of  uSr(x,t)  and 
V£r  (x)  =  s~2 111  w£r  (x,s)  E  C2+a  (Q)  be  the  corresponding  tail  function  (see  (15)).  Suppose 
that  the  cut-off  pseudo  frequency  s  is  so  large  that  for  any  such  function  er  (x)  the  following 
estimates  hold 

I  1 2+a  —  I^rl2 +a  (28) 

where  (  E  (0, 1)  is  a  sufficiently  small  number.  Denote 


p  \=  2  {h  +  cr  +  x  +  £) . 


(29) 


Let  K  be  the  constant  of  the  Schauder  theorem  in  (27)  and  N  <  N  be  the  total  number 
of  functions  qn  calculated  by  the  algorithm  of  subsection  3.1.  Suppose  that  the  number  N  = 
N  (h)  is  connected  with  the  step  size  h  via  N  (h)  h  =  f3,  where  the  constant  (3  >  0  is 
independent  on  h.  Let  f3  be  so  small  that 


13  ~  384/1  s2 


1 

24  KM*' 


(30) 


In  addition,  let  the  number  p  and  the  parameter  p  of  the  CWF  satisfy  the  following  estimates 


V<Vo  {K,C*,s) 


1 

16  KM* 


1 

256/1  C*s2  ’ 


(31) 


P>  Po  (■ C*,K,s,p ) 


48/1  C*s2, 
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Then  for  each  n  G  [l,lV]  the  sequence  {^1}^  converges  in  C2+a  (n)  to  the  function 
qn,\.  Likewise,  if  iterations  of  {qn,i}  with  respect  to  the  tails  are  stopped  at  i  =  mn  with 
q.n,mn  '■=  Qn  for  each  n  G  [l,iV]  ,  then  the  following  convergence  estimates  hold 

I  Qn  ~  Qn  1 2+a  —  M  ^ 

\qn\2+a  <  2 c*,  n  G  [l ,N]  , 

+  (33) 

3.3  Discussion  of  Theorem  1 

3.3.1  The  parameter  £ 

By  (16)  and  (28)  the  parameter  £  is  small  as  long  as  the  truncated  pseudo  frequency  s 
is  large.  This  implies  of  course  that  the  parameter  rj  in  (29)  is  also  small.  Nevertheless, 
this  theorem  has  a  discrepancy  related  to  the  parameter  £.  Indeed,  by  (31)  we  should  have 
V  <  Ci/s 2,  Ci  =  (256 KC*)-1 .  On  the  other  hand,  (16)  implies  that  £  =  O  (1/s) ,  s  — >  oo. 
Clearly  these  two  conditions  imposed  on  g  and  £  are  incompatible  with  (29).  In  addition 
since  by  (26)  M*  =  O  (s2)  as  s  — >  oo ,  then  there  is  no  guarantee  that  the  right  hand  side  of 
(32)  is  small  indeed. 

We  explain  this  discrepancy  as  follows.  Since  the  problem  of  construction  of  globally 
convergent  numerical  methods  for  our  CIP  is  obviously  an  extremely  challenging  one,  we 
need  to  make  certain  compromises  such  as,  e.g.  ones  outlined  in  the  previous  paragraph.  In 
simple  terms,  not  everything  can  be  covered  by  the  theory.  The  only  way  to  justify  these 
compromises  is  via  numerical  experiments.  Numerical  experiments  of  previous  publications 
on  this  method  [2,  3,  4,  5]  fully  confirm  the  theory,  see,  e.g.  subsection  9.3  in  [3].  Thus,  they 
demonstrate  that  this  compromise  is  reasonable.  Note  that  we  truncate  our  pseudo  frequency 
“gently”.  In  other  words,  instead  of  just  setting  for  the  tail  function  to  be  zero,  we  iterate 
with  respect  to  tails.  In  addition,  we  refer  to  the  Gelfand-Krein-Levitan  method  for  a  2-D 
inverse  hyperbolic  problem  of  [14].  This  method  shows  an  excellent  numerical  performance. 
Still,  it  has  a  similar  problem  with  the  truncation  of  the  Fourier  series  with  respect  to  one 
of  spatial  variables.  In  connection  with  this  compromise,  we  present  a  new  mathematical 
model  in  the  next  sub-subsection  and  we  also  demonstrate  that  a  very  similar  problem  takes 
place  in  the  classic  Real  Analysis.  Finally,  we  believe  that  results  of  the  current  publication 
eliminate  last  remaining  doubts  about  the  validity  of  the  technique  of  [2,  3],  see  section  8. 

3.3.2  A  new  mathematical  model  and  a  classic  example  from  Real  Analysis 

We  associate  with  Theorem  1  a  new  mathematical  model.  In  this  model,  as  soon  as  the 
large  truncation  pseudo  frequency  s  is  chosen,  we  allow  the  parameter  £,  which  bounds  tails 
in  (28)  and  which  is  involved  in  the  convergence  parameter  r)  in  (29),  to  be  infinitely  small, 
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independently  on  s.  Actually,  we  do  exactly  this  in  our  numerical  implementation,  which 
justifies  this  new  model,  see  subsections  7.1  and  7.2  for  details. 

Finally,  we  provide  here  an  example  of  the  same  nature.  This  example  is  linked  with  the 
classic  issue  of  Asymptotic  Series.  Although  these  series  quite  often  formally  diverge,  still 
their  truncations  provide  very  good  approximations  for  corresponding  functions.  Consider 
the  classic  error  function  4>  0*0  [!]> 


OO 


X 


The  asymptotic  series  expansion  for  the  function  ex2&  ( x )  is 


OO 


i  +  2.H0r 

n=l 


(2ra-  l)!l\ 

w ) 


,x 


OO, 


where  (2 n  —  1)!!  =  T3-5...-(2n  —  1) .  This  asymptotic  series  diverges  for  any  x.  Nevertheless, 
it  is  well  known  that  the  following  is  a  good  approximation  for  the  function  ex  $  (x)  for  large 
values  of  x, 


N 


(x) 


7 TX 


Ef-1)’ 


72—1 


(2n  —  1)!! 
(2x2)n 


,x 


OO. 


(34) 


The  truncation  number  N  in  (34)  has  exactly  the  same  nature  as  our  truncation  pseudo 
frequency  s.  Hence,  this  example  demonstrates  that  things,  which  are  similar  to  ours,  exist 
in  the  classic  Real  Analysis,  and  they  are  related  to  the  asymptotic  series  expansions. 


3.3.3  The  meaning  of  estimates  (30)  and  (33) 


The  estimate  (33)  tells  one  that  the  accuracy  of  the  reconstruction  of  the  function  e*  is 
improving  with  iterations  for  the  first  few  iterations.  However,  when  n  becomes  sufficiently 


large,  the  norm 


-H 


becomes  comparable  with  the  level  of  error  (23/8)  rj.  This  error 


includes  the  error  £  in  our  new  mathematical  model,  the  error  a  in  the  data  and  some  less 
significant  errors  h  and  x  of  our  method.  In  other  words,  we  cannot  guarantee  that  our 
accuracy  will  improve  after  reaching  a  certain  n  :=  no  G  [2,  TV]  .  In  addition,  by  (30)  we  also 
cannot  guarantee  anything  for  n  >  N .  We  indeed  observe  this  in  our  computations,  in  which 
we  take  either  N  =  6  or  N  =  5,  and  we  do  this  for  a  very  plausible  reason,  see  (45)  and 
subsection  7.2  below.  It  was  pointed  out  on  pages  156  and  157  of  the  book  [9]  that  one  of 
backbone  ideas  of  the  theory  of  Ill- Posed  Problems  is  to  use  the  number  of  iterations  as  one 
of  regularization  parameters.  This  iteration  number  is  N  in  our  case.  The  true  reason  why 
the  number  (5  is  small  in  (30)  is  that  equations  (17)  are  generated  by  equations  (13),  which 
contain  nonlinear  terms  with  Volterra  integrals.  It  is  well  known  that  one  can  guarantee 
existence  of  solution  of  a  nonlinear  Volterra  integral  equation  only  on  a  small  interval.  For 
example,  the  Cauchy  problem  y'  =  y2  +  l,y  (0)  =  0  has  its  solution  y  (z)  =  tan  z  with 
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the  singularity  at  z  =  tt/2.  On  the  other  hand,  this  Cauchy  problem  is  equivalent  with  the 
following  nonlinear  integral  equation  of  the  Volterra  type 

Z 

y(z)  =  I  y 2  (t)  dr  +  z. 

0 

4  Experimental  Setup 

4.1  Data  acquisition 


Source 


Figure  1.  Schematic  diagram  of  the  source/ detectors  configuration,  a)  The  prism  depicts  our  computational 
domain  f l.  This  domain  is  a  part  of  another  prism,  which  was  our  holder  made  out  of  Styrofoam.  Only 
a  single  source  location  was  used.  Tomographic  measurements  of  the  scattered  time  resolved  EM  wave  were 
conducted  on  the  bottom  side  of  this  prism,  b)  Schematic  diagram  of  locations  of  detectors  (probes)  on  the 
bottom  side  of  the  prism  f l.  The  distance  between  neighboring  probes  was  10  mm. 

Everywhere  below  the  word  “prism”  means  “rectangular  prism”,  a  conventional  notion 
of  geometry.  For  brevity  below  x  denotes  both  a  vector  x  G  M3  and  one  of  components 
of  this  vector  x  =  (x,  y,  z )  .  ft  is  always  clear  from  the  context  what  is  what  there.  Our 
source/detectors  configuration  is  schematically  depicted  on  Figure  1.  The  source  has  gener¬ 
ated  an  electromagnetic  (EM)  wave,  which  we  wanted  to  be  a  plane  wave  when  reaching  the 
bottom  side  of  the  prism  of  Figure  1,  where  measurements  were  conducted.  But  actually  this 
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was  a  spherical  wave,  because  of  a  rather  small  distance  between  the  source  and  that  side. 
We  had  a  holder  consisting  of  Styrofoam.  Styrofoam  is  a  material,  whose  relative  permittiv¬ 
ity  er  &  1,  i.e.  the  same  as  one  in  the  air.  Sizes  of  the  holder  were  260  mm  x  135  mm  x  260 
mm.  However,  because  of  our  previous  computational  experience  [3],  we  have  chosen  another 
prism  as  our  computational  domain  hi,  which  mostly  a  part  of  this  holder,  except  that  its 
size  in  the  y-direction  was  5  mm  more.  It  should  be  kept  in  mind  that  the  above  holder 
physically  existed,  whereas  £2  was  sort  of  “imaginary”  domain  which  we  have  used  for  our 
computations.  So  sizes  of  the  prism  hi  were  240  mm  x  140  mm  x  240  mm.  Hence,  sizes  of 
front  and  back  sides  of  the  prism  of  Figure  1  are  240  mm  x  240  mm,  sizes  of  other  four  sides 
are  240  mm  x  140  mm,  and  this  prism  is  exactly  our  domain  £2  in  (4).  The  distance  between 
the  wave  source  and  the  top  side  of  the  domain  £2  was  130  mm.  The  initializing  pulse  was 
100  ps  duration.  Since  the  speed  of  the  EM  wave  propagation  in  the  air  is  0.3  nnn/ps,  then 
it  requires  433  ps~130/03  ps  for  this  wave  to  travel  form  the  source  to  the  top  boundary  of 
£2.  Hence,  the  wave  did  not  yet  reach  the  prism  £2  during  the  100  ps  duration  of  this  pulse. 
The  initializing  pulse  was 


m 


A  sin  (^r)  ,  for  r  E  (0, 100)  ps, 
0,  for  r  >  100  ps, 


(35) 


where  A  is  the  amplitude.  Note  that  it  is  unclear  a  priori  on  how  the  value  of  A  will  be 
“reflected”  in  our  mathematical  model.  Still,  our  data  processing  procedure  does  not  rely 
on  a  knowledge  of  A.  The  time  resolved  signal  was  measured  at  some  locations  of  the  probe 
(i.e.,  detector)  on  the  bottom  side  of  the  prism  £2,  as  indicated  on  Figure  1-b). 

A  special  question  to  address  was  about  the  timing  of  measurements.  The  original  pulse 
needs  a  few  nanoseconds  to  go  through  the  device  and  reach  the  tip  of  the  EM  wave  generator, 
which  is  depicted  on  Figure  1  as  the  source  location.  In  our  mathematical  model  the  zero 
time  is  the  moment  when  the  pulse  leaves  the  tip  of  the  generator.  On  the  other  hand, 
on  each  location  of  the  probe  measurements  were  conducted  from  the  moment  when  the 
pulse  was  initiated  “within”  the  device,  which  is  prior  the  moment  when  it  leaves  the  tip 
of  the  generator.  Hence,  actual  measurements  were  conducted  for  times  r  E  (0, 12300)  ps 
=  (0, 12.3)  ns,  where  r  is  the  real  time  with  dimensions.  We  had  two  measurements  at  each 
probe  location:  (1)  First  we  have  measured  the  reference  signal  when  the  inclusion  was  not 
present,  and  (2)  Second,  we  have  measured  the  signal  when  the  inclusion  was  present. 

The  step  size  in  time  between  two  consecutive  measurements  was  A r  =  20  ps.  Hence, 
we  had  only  Eve  (5)  measurement  points  per  100  ps  duration  of  the  initializing  pulse.  We 
have  measured  the  scattering  EM  wave  sequentially.  In  other  words,  first,  we  were  putting 
the  probe  at  one  location,  sent  the  pulse  and  measured  the  time  resolved  scattering  wave 
at  this  location.  Next,  we  moved  the  probe  mechanically  in  a  neighboring  location  and 
repeated  the  measurement,  etc..  The  distance  between  two  neighboring  probe  locations 
was  10  mm,  and  so  we  have  covered  the  entire  bottom  side  of  the  holder  by  this  grid 
with  10  mm  step  size.  Pulses  were  generated  by  the  Picosecond  Pulse  Generator  10070A, 
see  http://www.picosecond.com/.  The  scattered  EM  wave  was  measured  by  Tektronix 
DSA70000  Series  Real  Time  Oscilloscope,  see  Tektronix,  http:/ /www. tek.com/products/oscilloscopes/. 
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a) 


b) 


Figure  2.  a)  Picosecond  Pulse  Generator  10070 A  and  b)  Tektronix  DSA70000  Series  Real  Time  Oscilloscope 


4.2  Dimensionless  variables 

To  work  with  the  data,  we  have  re-scaled  our  dimensions  in  time  and  space  and  have  made 
them  dimensionless.  First,  we  have  re-scaled  spatial  variables.  While  working  previously 
with  the  above  CIP  for  computationally  simulated  data  [3]  when  the  domain  fl  was  a  prism 
with  ratios  of  lengths  of  sides  similar  with  ones  of  Figure  1,  we  have  discovered  that  the  step 
size  h  =  0.2  when  solving  both  forward  and  inverse  problems  was  an  optimal  one.  On  the 
other  hand,  since  we  had  the  10  mm  distance  between  neighboring  positions  of  the  probe, 
we  have  decided  to  re-scale  spatial  variables  in  such  a  way  that  10  mm  would  turn  into  0.2  in 
dimensionless  variables.  So,  let  ( x',y',z ')  be  spatial  variables  with  dimensions.  Then,  since 
10/02=50,  we  set  for  re-scaled  dimensionless  variables  (. x,y,z )  =  (x,,y,,zl)  / 50.  Thus,  our 
dimensionless  computational  domain  fl  for  the  CIP  is 

«  =  {{x,y,z)  e  [-2.4, 2.4]  x  [-1.4, -1.4]  x  [-2.4, 2.4]}.  (36) 

Denote  P  the  bottom  side  of  the  domain  fl  in  (36), 

P  =  {(x,  y,  z)  :  (x,  y)  G  [-2.4,  2.4]  x  [-1.4,  -1.4],  z  =  -2.4}  .  (37) 

We  now  should  address  the  question  on  how  long  we  should  measure  the  scattered  EM 
wave  on  probes  located  on  P.  It  was  not  easy  to  find  out  in  our  experimental  arrangement 
when  exactly  the  pulse  leaves  the  tip  of  the  EM  waves  generator,  i.e.  the  source  depicted  on 
Figure  1.  However,  we  knew  that  the  signal  arrives  at  the  probe  approximately  at  11,520 
ps.  Since  the  distance  between  the  planar  surface  P  and  the  source  was  370  mm,  the  speed 
of  light  in  the  air  is  0.3  mm/ps  and  (370mm)  /  (0.3 mm/ps)  =  1233  ps,  then  the  zero  time 
should  be  at  11,  520  ps— 1,  233  ps  ~  10,  300  ps:=  r0.  So,  we  should  work  with  a  new  variable 
t'  =  t  —  Tq.  The  next  question  was  on  how  to  re-scale  the  time  r'.  It  follows  from  (1),  (3)  and 
(4)  that  the  refractive  index  outside  of  the  domain  D  is  n  (x)  =  1.  This  means  that  the  EM 
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wave  should  travel  the  dimensionless  distance  of  h  =  0.2  in  0.2  dimensionless  time  units.  On 
the  other  hand,  0.2  corresponds  to  10  mm.  Let  t  denotes  the  dimensionless  time.  Then  we 
should  choose  such  a  multiplier  7  >  0  (7  has  dimension  in  picoseconds)  that  7 1  =  r' .  Hence, 
we  should  have  O.27 ps  =  10 mm/  (0.3 mm/ps) ,  which  implies  that  7  =  166.67  ps.  Thus, 

7 1  =  t' ,  7  =  166.67  ps.  (38) 

Finally,  we  should  figure  out  on  how  long  we  should  measure  the  output  signal.  In  [3]  we 
have  worked  with  the  time  interval  t  G  (0,T)  =  (0, 12) .  Since  by  (38)  12y  fa  2000  ps,  then 
we  should  work  with  t'  G  (0,  2000)  ps.  So,  since  by  the  above  r0  =  10,  300  ps  and  r'  =  r  —  r0, 
then  the  maximal  value  of  r  of  our  interest  is  12,  300  ps  =  10,  300  ps  +2000  ps.  Hence,  we 
should  measure  the  output  signal  for  r  G  (0, 12300)  ps  =  (0, 12.3)  ns. 

Remark.  When  re-scaling  in  this  sub-subsection,  we  have  not  reflected  this  in  the  cor¬ 
responding  PDE.  This  becomes  possible  because  of  the  data  processing  procedure  described 
in  section  6.  Below  we  work  only  with  dimensionless  time  and  spatial  variables. 

4.3  Why  measuring  the  reference  signal  at  each  probe  location 
rather  than  at  a  single  one 

In  principle,  our  technique,  including  the  data  processing  described  below,  allows  the  mea¬ 
surement  of  the  reference  signal  only  at  one  probe  location  outside  of  the  medium  of  interest: 
for  calibration  purposes.  The  only  reason  why  we  have  measured  the  reference  signal  for  each 
location  of  the  probe  was  that  our  current  numerical  implementation  of  the  globally  con¬ 
vergent  algorithm  works  only  with  the  case  when  the  initializing  wave  field  is  a  plane  wave. 
On  the  other  hand,  a  visual  inspection  of  the  output  experimental  data  has  revealed  to  us 
that  actually  we  had  a  spherical  rather  than  a  plane  wave.  An  extension  of  our  numerical 
implementation  to  the  case  of  the  spherical  wave  is  rather  straightforward  and  we  plan  it 
for  the  future.  We  have  used  the  point  source  rather  than  a  plane  wave  in  (1),  (2)  only 
to  obtain  the  asymptotic  behavior  (16),  which  actually  follows  from  the  construction  of  the 
fundamental  solution  of  the  hyperbolic  equation  in  [20,  21].  As  it  was  mentioned  in  section 
2,  in  computational  practice  we  verify  this  asymptotic  behavior  numerically  when  working 
with  plane  waves,  see  subsection  7.2  of  [2], 

5  Data  Processing 

5.1  Data  simulation 

Since  the  computationally  simulated  data  play  an  important  role  in  our  data  processing 
procedure,  we  first  describe  the  solution  of  the  forward  problem  for  equation  (1).  We  have 
computationally  simulated  the  solution  of  the  forward  problem  for  the  reference  medium. 
Since  it  is  impossible  to  actually  solve  the  PDE  (1)  in  the  entire  space  M3,  we  have  solved 
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it  in  a  such  a  domain  G  that  O  C  G,  where  G2  is  defined  in  (36).  We  took  the  domain 
G  =  {(x,  y1  z)  G  [—3,  3]  x  [—2,  2]  x  [—5, 5]}  .  Our  initializing  plane  wave  was  v  ( t ) , 

[  sin  (ut)  ,  for  t  G  (0,  A)  , 
v(t)  =  l  0,  for  t  >  A , 

[  u  =  7. 

Let  dG i  and  8G2  be  respectively  top  and  bottom  sides  of  the  prism  G  and  dG3  =  dG\  (dG i  U  dG2 ) 
be  the  rest  of  the  boundary  of  the  domain  G.  We  have  numerically  solved  the  following  initial 
boundary  value  problem 


er  (x)  Utt  =  A  u,  in  Gx  (0  ,T),T  =  12, 

u(x,  0)  =  0,  ut(x,  0)  =  0,  in  G, 

,  27 r 

dnu\aCl  =  v  (t) ,  on  dG  1  x  (0,  — 

dnu\dGi  =  ~dtu,  on  dGi  x  (h,T), 

dnu\ ac2  =  -Ot'u.  on  dG2  x  (0 ,T), 

dnu\ 9G3  =  0,  on  dG3  x  (0,T), 


(39) 


In  the  case  of  data  simulation  for  the  reference  medium  we  have  in  (39)  er  (x)  =  1.  We  denote 
this  solution  as  uref  (x,  t ) .  Thus,  in  (39)  the  plane  wave  is  initialized  at  the  top  boundary 
3G\  for  t  G  (0,  27t/7]  and  propagates  into  G  .  First  order  absorbing  boundary  conditions 
[10]  were  used  on  the  top  boundary  for  t  G  (27t/7,  T)  and  on  the  bottom  boundary  dG2  for 
t  G  (0,  T).  The  zero  Neumann  boundary  condition  is  used  on  the  rest  of  the  boundary  of  the 
prism  G.  The  latter  boundary  condition  is  used  because  the  “pure”  plane  wave  with  er  (x) 
=  1  satisfies  this  condition.  The  problem  (39)  was  solved  by  the  hybrid  FEM/FDM  method 
described  in  [6].  In  this  method,  FDM  is  used  outside  of  the  domain  fl  and  FEM  is  used 
inside  of  this  domain.  The  step  size  in  the  overlapping  region  was  h  =  0.2,  i.e.  this  was  the 
same  step  size  as  the  distance  between  neighboring  probes. 


5.2  Measured  time  resolved  data 

Denote  xm  G  P  the  location  of  the  probe  number  m  at  the  bottom  side  P  of  the  prism 
fl,  see  (36)  for  P.  Figures  3-a)-d)  display  samples  of  curves  and  the  label  for  Figure  3 
explains  details.  The  first  arrival  signal  is  when  the  burst  starts.  The  signal  before  it  is  a 
low  frequency  signal  reflecting  some  processes  in  the  probe  itself.  It  is  obvious  that  the  first 
thing  to  do  is  to  perform  the  Fourier  transform  and  truncate  too  low  and  too  high  frequencies. 
Low  frequencies  should  be  truncated  to  diminish  the  signal  resulting  from  the  probe  itself. 
And  high  frequencies  should  be  truncated  to  somehow  decrease  the  noise  further.  We  are 
interested  in  the  first  burst  only.  We  see  on  Figs.  3-c),d)  that  the  signal  before  this  burst 
looks  like  a  statistical  noise.  Even  if  it  is  not  exactly  a  statistical  noise,  we  have  decided 
to  subtract  it  from  the  data  in  the  time  interval  where  the  first  burst  is.  Let  to  be  the 
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Figure  3.  a)  A  sample  of  the  measured  reference  time  resolved  signal  (i.e.,  no  inclusion  present)  at  the 
location  xm  £  P  of  the  probe  number  m.  b)  The  measured  signal  with  inclusion  present  at  the  same  probe 
location.  The  first  burst  starts  when  the  EM  wave  arrives  at  the  probe.  The  signal  before  this  burst  reflects  a 
process  within  the  probe  itself,  c)  and  d)  represent  signals  a)  and  b)  respectively  after  cleaning  some  noise  via 
applying  the  Fast  Fourier  Transform  procedure  of  MATLAB  and  truncating  too  low  and  too  high  frequencies. 
We  are  interested  in  the  area  of  the  first  burst  only.  One  can  observe  that  the  amplitude  of  the  signal  with 
the  dielectric  inclusion  present  (Fig.  3d)  is  generally  less  than  one  of  the  reference  signal. 

approximate  time  when  the  first  burst  starts  and  t\  be  the  approximate  time  when  this 
burst  ends.  Consider  an  interval  (a,  b )  with  0  <  a  <  b  <  t0  such  that  b  —  a  —  t\  —  to-  For 
either  of  Figures  3-c)  or  3-d)  let  fa,b  (t)  be  the  signal  on  this  interval,  and  let  the  first  burst 
be  described  by  the  function  f  (t)  ,t  G  (t0,ti).  We  have  replaced  f  (t)  with  the  function 
/  (t)  =  f  (t)  —  fa,b  (t  —  t o  +  a)  for  t  G  (to,  ti)  and  have  worked  with  this  function  only  after 
this.  Figure  4  displays  of  resulting  superimposed  curves  for  both  reference  signal  and  the 
signal  with  inclusion  present.  These  curves  are  generated  by  Figures  3-c),  d).  We  have  also 
set  to  zero  those  parts  of  these  curves  which  were  before  the  first  burst. 
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Amplitude 


Super  Imposed  of  Reference  Measurement  and  Inclusion  Measurement  Super  Imposed  of  Representation  with  Model  Curve  of  Reference  Measurement  and  Inclusion  Measurement 


a) 


b) 


Figure  4.  This  figure  explains  the  idea  of  the  immersing  procedure  in  the  time  domain,  a)  Resulting  superim¬ 
posed  experimental  curves  obtained  from  curves  on  Figures  3-a),  b).  The  red  curve  is  for  the  reference  signal 
and  the  blue  curve  is  for  the  signal  with  a  dielectric  inclusion  present,  both  at  the  same  location  xm  £  P 
of  the  probe  number  m.  b)  The  red  curve  displays  computationally  simulated  data  uref  (xm,t),  also  see 
sub-subsection  6.3.2.  The  blue  curve  Uinci  ( xm,t )  =  uref  ( xm,t  —  A tm)  K™p/M™p  represents  a  sample  of  the 
immersed  experimental  data  in  the  time  domain  at  the  same  probe  location  xm  £  P,  see  explanations  in 
sub-subsection  6.f.l.  It  is  only  the  blue  curve  with  which  we  work  further.  The  red  curve  is  displayed  for  the 
illustration  purpose  only. 


5.3  The  key  procedure:  immersing  in  the  computationally  simu¬ 
lated  data 

Consider  now  Figure  4-a)  for  a  probe  xm  G  P  number  m.  We  have  decided  to  “immerse”  our 
experimental  data  in  the  computationally  simulated  data  using  the  largest  peak  in  the  red 
curve  (reference  medium)  with  the  peak  value  M™p  and  the  next  peak  after  it  in  the  blue 
curve  (the  medium  with  a  dielectric  inclusion  present)  with  the  peak  value  Kf)p-  This  next 
peak  was  chosen  because  the  presence  of  a  dielectric  inclusion  results  in  a  time  delay  of  the 
EM  wave.  This  idea  led  us  to  a  radically  new  data  processing  procedure  described  in  this 
subsection.  The  immersing  procedure  in  the  computationally  simulated  data  consists  of  two 
stages:  immersing  in  the  time  domain  and  subsequent  immersing  in  the  pseudo  frequency 
domain.  We  describe  both  of  them  sequentially.  Below  we  talk  about  the  computational 
simulated  data  in  several  places.  In  all  cases  these  simulated  data  were  computed  before  the 
experimental  data  were  measured. 

5.3.1  Stage  1.  Immersing  in  time  domain 

Recall  that  the  function  uref  ( x ,  t )  is  the  solution  of  the  problem  (39)  with  computationally 
simulated  data  for  er  =  1.  Obviously  uref  (x^ft)  =  uref  ,  Va^1),  x^  G  P,Vt  G  (0,T) . 

Let  t  :=  tsr g™  be  such  a  moment  of  time  that  for  all  x  G  P  we  have  uref  ( x ,  t)  =  0  for  t  <  tffj 
and  uref  ( x ,  t)  >  0  for  such  moments  of  time  t  that  are  rather  close  to  tsffff  with  t  >  tffff,  see 
the  reference  curve  (red)  on  Figure  4-b).  Naturally,  we  call  t\ fffff  the  time  of  the  first  arrival 
of  the  computationally  simulated  EM  wave  field  u  ( x ,  t )  =  uref  ( x ,  t )  in  (39)  at  the  plane  P. 

Important  Observation  and  Conclusion.  We  have  observed  in  our  computational 
simulations  that  for  those  values  of  t  >  tsffj,  which  are  rather  close  to  tsrlff,  the  only  side  of 
which  is  sensitive  to  the  presence  of  a  dielectric  inclusion,  is  the  bottom  side  P  of  the 
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domain  f2.  Other  five  sides  of  the  prism  0  are  not  sensitive  to  the  presence  of  inclusions  for 
those  values  of  t.  In  other  words,  the  values  of  the  function  u  at  those  sides  are  approximately 
the  same  ones  as  those  for  the  reference  medium.  This  important  observation  led  us  to  the 
conclusion  that  we  should  work  only  with  the  information  contained  in  the  first  burst,  see 
Figure  4-a).  Actually,  this  idea  corresponds  well  with  the  common  knowledge  of  physics  that 
the  most  informative  signal  is  the  one  which  is  collected  at  those  times  t,  which  are  close  to 
the  moment  of  the  first  arrival  of  the  signal. 

The  next  important  question  was:  How  to  work  with  this  first  burst?  Indeed,  we  have 
observed  on  experimental  curves  that  the  amplitude  of  the  signal  near  the  time  of  the  first 
arrival  at  a  probe  is  so  weak  that  it  this  signal  cannot  be  differentiated  from  the  noise,  see 
Figure  4-a).  Consider  the  reference  experimental  signal  at  the  probe  xm  G  P,  see  Figure  4-a). 
We  have  decided  that  we  should  work  with  the  largest  upwards  looking  peak.  If,  however  we 
have  several  peak  values  with  no  more  than  10%  difference  between  them,  then  we  choose 
the  earliest  among  those,  i.e.  we  choose  such  a  peak  which  corresponds  to  the  minimal 
value  of  t.  We  have  done  this  as  follows.  Let  ym  =  ym  (t)  be  the  experimentally  measured 
function  at  the  point  xm  (after  the  above  noise  reduction).  Suppose  that  on  the  first  burst 
we  have  local  maxima  at  points  {H, ...,  tn}  and  their  values  are  y±  =  ym  {tfi) , ..,  yn  =  ym  ( tn ) . 
Let  Y  =  ym  (tk)  =  maxi=i (U)  ■  Consider  numbers  y1  =  y  (t1) ,  ...,yr  =  y(tr ) ,  where 
points  {f1,...,^}  C  {ti,...,tn}  are  such  that  y3  /  Y  G  [0.9,1].  And  let  tm  =  min  {t1, ...,  tr}  . 
We  choose  the  pair  ( tm,ym )  and  denote  ( tm,ym )  :=  ( m ) ,  M™p)  . 

Next,  we  choose  a  local  maximum  for  the  experimental  curve  at  {xm}  for  the  medium 
with  a  dielectric  inclusion  present.  Let  zm  =  zm  ( t )  be  that  curve  (after  the  above  noise 
reduction) .  Consider  local  maxima  of  this  function  for  t  >  ( m ) .  Let  tefifificl  ( m )  >  (m) 

be  the  minimal  value  of  the  time  t  on  the  interval  {t  >  (m)}  at  which  a  local  maximum 

of  the  function  zm  ( t )  is  achieved.  In  other  words,  we  choose  the  first  upwards  looking  peak 
of  the  function  zm  ( t )  occurring  after  the  prior  chosen  reference  peak  at  the  reference  curve. 
Denote  zm  M)  :=  K V  So,  KZr  is  the  value  of  that  peak,  see  Figure  4-a).  However,  if 
^exp/Mexp  >  2  A  then  we  set  for  the  point  {xm}  that  (t^,  (m) ,  K™xp)  :=  (m) ,  M™p)  . 

We  have  observed  that  on  all  probes  /T™p  <  M™p. 

Now  we  are  ready  to  immerse  our  data  in  the  computationally  simulated  data.  Let  u  ( x ,  t) 
be  the  solution  of  the  problem  (39)  with  er  ( x )  =  1.  Hence,  the  function  u  ( x ,  t )  describes  the 
plane  wave  propagation  in  the  uniform  medium.  Let  xm  G  P  be  the  location  of  the  probe 
number  m.  Let  A tm  =  tefifificl  ( m )  —  (m)  be  the  time  shift  between  two  above  chosen  peaks 

for  this  probe,  see  Figure  4-a).  Then  we  set 

Km 

Uincl  (xm,  t )  =  -TJ^-Uref  (xm,  t  ~  A tm)  .  (40) 

2^exp 

So,  (40)  is  our  immersed  data  in  the  time  domain  for  the  probe  number  m.  Figure  4-b) 
illustrates  (40).  It  is  clear  from  the  above  that  if  K™p/M™p  >  2/3,  then  uinci  (xm,t)  = 
uref  (xm,t) .  Below  we  work  only  with  so  immersed  data  for  each  probe  location  on  the 
surface  P.  We  work  further  with  these  immersed  data  to  ultimately  use  them  as  an  input  for 
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the  Dirichlet  boundary  conditions  i\)n  (x)  ,x  £  P  in  equations  (17).  Thus,  since  we  actually 
use  ratios  iF™p/M™p,  we  do  not  need  to  know  the  value  A  of  the  signal’s  amplitude  in  (35). 


5.3.2  Stage  2.  Immersing  in  the  pseudo  frequency  domain 


We  should  apply  the  Laplace  transform  (7)  to  each  curve  uinci  (. xm ,  t)  in  (40)  to  obtain  the 
function  Wj.nci  (xm,  s)  =  C  ( Uinci  (xm,  t)).  The  next  question  is:  For  what  values  of  the  pseudo 
frequency  s  should  we  actually  calculate  the  integral  (7)?  To  address  this  question,  we  have 
solved  the  above  CIP  for  the  domain  f2  in  (36)  for  a  computationally  simulated  data  with 
an  inclusion  present,  using  the  above  globally  convergent  numerical  method.  This  was  done 
prior  obtaining  the  experimental  data.  We  have  established  that  the  following  numbers  are 
optimal  ones 

s  €  [3.5,  7.5] ,  h  =  0.5.  (41) 

So,  we  have  calculated  the  numbers  Winci  (. xm ,  s )  for  nine  values  of  s  in  (41)  for  each  probe 
location  xm  €  P.  Below  in  this  sub-subsection  we  work  only  with  values  of  s  from  (41).  We 
have  observed  computationally  that,  because  of  the  rapid  decay  of  the  function  exp  (— st ) 
with  respect  to  t  for  values  of  s  from  (41),  the  major  impact  in  the  integral  (7)  comes  from 
the  first  splash  of  the  curve  uref  (xm,t) .  So,  only  this  splash  is  depicted  on  Figure  4-b). 

Let  Wind(xm,s )  =  —  (In  Winci(xm,s))  / s2  and  for  each  value  s  from  (41)  let  Winci(x,s ) 
be  the  linear  interpolation  of  discrete  values  {winci  (xm,  s)}  over  P.  The  function  winci  (x,  s ) 
is  very  noisy  with  respect  to  x  £  P,  e.g.  see  Figure  5- a).  On  the  other  hand,  Figure  5-b) 
displays  a  typical  x—  dependence  of  the  function  wsim  (x,s)  :=  —  (In  wsim  (x,s))  /s2,  x  £  P, 
where  the  function  wSim  (x,  s )  is  corresponds  to  a  sample  of  computationally  simulated  data 
for  an  inclusion  present.  Again,  this  data  simulation  took  place  prior  the  experimental  data 
were  obtained.  Hence,  to  make  our  resulting  function  look  like  the  one  on  Figure  5-b),  we 
have  smoothed  the  function  winci  (x,  s )  over  x  £  P  via  the  so-called  Lowess  Fitting  in  the  2-D 
case,  which  we  took  from  MATLABR  2009a.  A  sample  of  the  resulting  function  wsmooth  (x,  s) 
for  s  :=  s  is  displayed  on  Figure  5-c).  Still,  however,  this  function  does  not  yet  look  similar 
to  the  function  depicted  on  Figure  5-c).  Define  the  number  Wref  (s)  as  follows  (see  (7)) 
Wref  (s)  =  — s~2  In  C  ( uref ) ,  for  x  £  P.  We  took 


w. 


immers 


(x,s) 


Wsmooth  (x,  s ) ,  if  iusmooth  (x,  s )  >  0.985  maxp  (■ wsmooth  (x,  s)) , 

Wref ,  otherwise, 


see  Figure  5-d),  which  depicts  the  function  Wimmers  {x,s) .  We  call  the  function  Wimmers  (x,  s) 
the  immersed  data  in  the  pseudo  frequency  domain.  Thus,  we  work  with  the  function 
Wimmers  (x,  s)  to  obtain  functions  ifn(x),x  £  P  in  (17).  Namely,  we  use  finite  differences 
to  approximately  compute  the  s-derivative.  Recalling  that  w  (x,  s )  :=  —  (lnw  (x,  s))  /s2,  we 
obtain  for  the  finite  difference 


x  = 


Wimmers  {x^  Sn  0.5)  Wimmers  (x j  SnJ 

(L5 


,  x  £  P. 


(42) 


As  to  the  values  of  the  function  ifn  (x)  on  other  five  sides  of  the  prism  f2,  they  were  computed 
by  the  same  finite  difference  formula  using  the  function  resulting  from  the  function  C  ( uref )  • 
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Showing  the  Last  Step  of  Post  Processing  at  s=7 


Showing  Fitting  Value  of  Laplace  Trasform  of  u(x,y,t)  on  Bottom  Face  at  s=7 


<5 


C)  d) 

Figure  5.  a)  The  function  Winci  (x,  s),  s  =  7.5.  b)  The  function  —  (In  wSim{x,s))  /s2  is  depicted,  where 
Wsim{x,s)  is  the  Laplace  transform  (7)  of  the  function  uSim{x,t)  for  a  computationally  simulated  data. 
Figure  5-b)  is  given  only  for  the  sake  of  comparison  with  Figure  5-a).  c)  The  function  wsrnooth  (x,  s)  resulting 
from  fitting  of  a)  by  the  Lowess  Fitting  procedure  in  the  2-D  case,  see  MATLABR  2009a.  d)  The  final  function 
Wimmers  (x,  s) .  Values  of  Wimmers  (x,  s)  are  used  to  produce  the  Dirichlet  boundary  conditions  i/jn  (x)  for 
PDEs  (17)  of  the  globally  convergent  algorithm,  see  (f7). 


5.3.3  Physics  considerations  for  data  immersing 

The  above  immersing  procedure  in  the  time  domain  makes  sense  from  the  physics  standpoint. 
Since  our  assumption  is  that  the  relative  permittivity  in  the  domain  of  interest  £7  is  er  >  1, 
then  by  (3)  the  EM  wave  should  propagate  through  slower  than  through  the  air.  In  other 
words,  the  so-called  time  delay  should  take  place,  and  this  is  why  (40)  makes  sense.  Still,  it 
is  unclear  how  to  actually  quantify  the  previous  statement.  Furthermore,  we  doubt  that  it  is 
known  what  exactly  that  statement  means  in  terms  of  the  experimental  data.  Does  it  relate 
to  the  truly  first  arrival  signal,  which  is  too  weak  to  be  trusted?  Or  perhaps  it  is  related  to 
the  largest  peak  values  as  we  have  chosen?  Or  maybe  this  is  related  to  the  first  burst  “as  a 
whole”? 

However,  if  we  assume  that  the  first  largest  peak  value  (although  within  10%  tolerance 
level)  for  the  reference  medium  arrives  earlier  for  the  air  than  the  one  for  the  case  of  a 
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heterogeneous  medium,  then  the  above  choice  of  peak  values  makes  sense.  Our  blind  imaging 
results  justify  the  latter  assumption.  Also,  using  only  one  point  ( tm ,  ym)  for  the  reference 
signal  and  only  one  point  zm  {t™d))  of  the  signal  with  an  inclusion  present  does  not 
mean  of  course  that  we  have  actually  used  only  this  single  point.  Indeed,  to  figure  out  which 
exactly  point  should  be  chosen,  we  have  counted  all  local  maxima  in  the  first  burst  and  thus, 
have  examined  the  entire  curve. 


6  Some  Details  of  the  Numerical  Implementation  of 
the  Globally  Convergent  Algorithm 

We  point  out  that  all  details  of  the  numerical  implementation  of  the  globally  convergent 
algorithm,  which  are  described  in  this  section,  were  implemented  prior  obtaining  the  exper¬ 
imental  data.  The  same  is  related  to  all  results  for  computationally  simulated  data  mentioned 
in  the  above  sections.  When  working  with  the  experimental  data,  we  have  not  changed  our 
original  numerical  code  for  the  inverse  problem  solution  and  thus  have  not  changed  features 
listed  in  the  next  subsection.  In  other  words,  our  computations  of  images  from  experimental 
data  were  unbiased.  We  have  implemented  all  these  details  listed  in  subsection  7.1  when 
working  with  those  computationally  simulated  data.  When  implementing  these  details,  our 
goal  was  twofold:  (1)  to  obtain  the  best  performance  of  the  globally  convergent  algorithm, 
and  (2)  at  the  same  time,  still  to  remain  “within”  conditions  of  Theorem  1. 


6.1  Details 


First  of  all,  it  is  important  to  point  out  that  all  details  listed  in  this  subsection  were  im¬ 
plemented  when  we  were  working  with  computationally  simulated  data,  which  was  prior 
obtaining  experimental  data. 

We  have  observed  in  our  computational  simulations  that  values  of  \Vn^  (x,s)|  of  approxi¬ 
mations  for  the  tail  function  dominate  values  of  all  other  terms  in  equations  (17)  for  functions 
qntk-  Hence,  when  solving  equations  (17)  for  functions  qnji  (subsection  3.1)  in  our  compu¬ 
tations,  we  have  used  in  (17)  s-derivatives  of  tails  djVn!i(x,s)  instead  of  tails  VU)i(x,s) 
themselves.  These  derivatives  were  taken  via  finite  differences,  similarly  with  (43).  How¬ 
ever,  when  computing  functions  (x)  via  (18a, b),  we  have  still  used  in  (18b)  the  function 
Vn,t  ( x ,  s)  itself.  The  parameter  in  the  Carleman  Weight  Function  was  fi  =  50.  Likewise,  we 
have  made  cut-offs  of  computed  functions  £r'  (x)  as  follows.  For  each  n  we  have  chosen 
a  cut-off  value  Ccut  ( n )  such  that  we  have  assigned  a  new  value  (x)  for  the  function 
Sr  (x)  as 


£rl,l)  (x) ,  if  e[n'l)  (x)  >  1  +  Ccut  (n) , 
1,  if  e^'l)  (x)  E  [1, 1  +  Ccut  (n)] . 


(43) 
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Note  that  by  (18a)  £r'’%'>  (x)  >  1,  Vx  E  Q.  The  numbers  Ccut  were  chosen  as  follows 

Ccut  (1)  =  0,  Ccut  (2)  =  0.2,  Ccut  (3)  =  C Cut  (4)  =  0.8,  Ccut  (5)  =  0.6, 

Ccut  (6)  =  Ccut  (7)  =  0.4,  Ccut  (8)  =  0.8. 

We  now  define  stopping  rules  of  iterations  for  functions  1  with  respect  to  the  nonlinear 
term  as  well  as  for  functions  {qn,k}  with  respect  to  the  tails.  These  rules  are  almost  the  same 
as  ones  in  [3].  Consider  the  planar  surface  P t  which  is  parallel  to  the  surface  P  in  (37).  The 
surface  P r  is  obtained  from  the  surface  P  via  shifting  upwards  by  h  =  0.2, 


P~h  =  {(; x,y,z)  :  (x,y)  E  [-2.4,  2.4]  x  [-1.4, -1.4],  «  =  -2.4  +  h  =  -2.2}  . 


And  let  IT  =  {(x,y)  E  [—2.4, 2.4]  x  [—1.4,  —1.4]}  be  the  orthogonal  projection  of  both  sur¬ 
faces  P  and  P^  on  the  (x,y)  plane.  Consider  norms  =  ||g^i|p-  —  ipn\\ l2(ci')-  We  stop 
iterations  of  functions  q ^  1  when  either  F%+1  >  F%,  or  norms  F%  stabilize,  or  F%  <  £,  where 
e  =  0.001  is  a  small  tolerance  number  of  our  choice.  Next,  we  iterate  with  respect  to  the 
tails.  We  similarly  introduce  norms  Fni  =  ||gn)j|^  —  V’nll L2(n')  an<4  use  the  same  stopping 
rule  as  one  for  F*. 

While  stopping  rules  for  iterations  of  functions  q *  1  and  qn,t  was  the  same  as  one  in  [3] , 
the  stopping  rule  for  computing  functions  £ ^  (x)  is  now  different  from  [3].  Namely,  let 
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bn  E  [1.9,4]  and  n  >  3,  then  take  the  final  solution  =  £^\ 
bn>  4  and  n  >  3,  then  take  the  final  solution  e£  = 
alternatively  compute  £r"+1^ 


(44) 


We  have  chosen  n  >  3  in  (44)  because  we  have  observed  in  our  work  with  computationally 
simulated  data  that  images  are  becoming  more  or  less  close  to  the  correct  ones  only  starting 
from  n  =  4.  And  this  is  in  a  conjunction  with  the  convergence  estimate  (33). 


6.2  Analytical  aspects  of  details  of  the  numerical  implementation 

Consider  the  analytical  aspect  of  details  listed  in  subsection  7.1.  Note  that  the  exact  mean¬ 
ing  of  words  “large”  and  “small”  (in  terms  of  used  numbers)  depends  on  a  specific  numerical 
implementation,  and  so,  from  this  standpoint,  we  can  assume  that  the  s-interval  (41)  corre¬ 
sponds  to  large  values  of  s  and  the  interval  (41)  is  small,  ft  follows  from  (16)  that  one  should 
expect  that 

\d-gV  (x,  s)|2+q,  <<  \V  (^,s)|2+q.  (45) 
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Denote  VU)l  (x,  s )  = 
of  the  right  inequality  (28), 


[x,  s )  .  It  follows  from  (45)  that  for  large  s  we  have  an  analog 


(*D  ’®) 


<  £  for  a  small  parameter  £,  whereas  the  left 

2+o 


inequality  (28)  is  still  valid.  It  easily  follows  from  the  proof  of  the  global  convergence 
theorem  in  [3]  that  Theorem  1  is  still  valid  when  tails  Vn/l  (x,  s)  are  replaced  with  Vrhl  (x,  s) . 
Therefore,  the  replacement  of  Vn^  (x,  s)  with  dgVnii  (x,  s)  corresponds  well  with  the  above  new 
mathematical  model  (sub-subsection  3.3.2).  This  is  because  by  this  model  one  considers  the 
number  £,  which  bounds  norms  of  tails  from  the  above,  as  a  parameter,  which  can  be  made 
infinitely  small  independently  on  the  value  of  s.  So,  we  conclude  that  in  our  computations  we 
have  indeed  used  the  above  new  mathematical  model.  As  to  the  cut-offs  (43),  such  cut-offs 
are  routinely  used  in  imaging  as  post  processing  procedures. 


7  Imaging  Results 

7.1  Abnormalities  and  their  positions 

Our  abnormalities  to  be  imaged  were  two  wooden  cubes,  see  Table  1. 


Table  1.  Sizes  of  two  wooden  cubes 


Cube  number 

Original  sizes,  mm 

Dimensionless  sizes 

1 

2 

40  x  40  x  40 

60  x  60  x  60 

0.8  x  0.8  x  0.8 

1.2  x  1.2  x  1.2 

Let  CL  be  the  center  line,  i.e.  the  straight  line  which  is  orthogonal  to  the  plane  P  and 
which  passes  through  the  source  of  EM  waves  (Figure  1).  Then  CL  —  {(x,y,  z)  :  x  =  y  =  0}  . 
We  have  placed  both  those  cubes  in  two  positions.  In  the  first  position  the  center  of  the 
cube  was  on  CL.  In  the  second  position  the  center  of  the  cube  was  shifted  off  CL  by  10  mm 
in  the  positive  direction  of  x-axis  (0.2  in  dimensionless  units).  In  addition,  we  have  used  the 
third  position  for  the  cube  number  1.  This  position  was  quite  off  CL.  Namely,  in  the  third 
position  the  center  of  this  cube  was  shifted  off  CL  by  60  mm  in  the  positive  direction  of 
the  rc-axis  (1.2  in  dimensionless  variables).  We  have  observed  on  the  experimental  data  that 
since  we  had  a  spherical  rather  than  a  plane  wave,  the  magnitude  of  the  EM  field  has  decayed 
quite  significantly  when  the  point  has  moved  rather  far  off  CL.  So,  the  goal  of  placing  the 
cube  number  1  in  the  third  position  was  to  see  how  this  decay  of  the  magnitude  of  the  EM 
held  would  affect  the  image  quality.  Because  of  some  logistical  reasons,  we  did  not  put  cube 
number  2  in  this  third  position.  Likewise,  because  of  logistical  reasons,  we  have  measured 
the  scattering  held  from  cube  number  1  in  the  hrst  position  twice:  in  two  consecutive  days. 
Therefore,  we  have  obtained  total  six  (6)  pieces  of  data  for  the  case  of  an  inclusion  present. 
In  addition,  the  data  for  the  reference  medium,  was  measured  only  once.  Table  2  lists  all  six 
cases. 
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Table  2.  Positions  of  centers  of  two  wooden  cubes  to  be  imaged  in  six  cases.  The  difference  between  cases 
1.1(1)  and  1-1(2)  is  that  they  were  measured  on  two  different  days  for  the  same  position  of  cube  1 


Cube  Number 

Case  Number 

Center 

1 

1.1(1) 

(0,0, -1.2) 

1 

1.1(2) 

(0,0, -1.2) 

1 

1.2 

(0.2,0, -1.2) 

1 

1.3 

(1.2,0, -1.2) 

2 

2.1 

(0,0, -1.2) 

2 

2.2 

(0.2,0, -1.2) 

7.2  Tables  and  images 

We  have  made  computations  using  the  above  described  globally  convergent  algorithm  and 
functions  (x)  in  (42).  We  have  used  the  stopping  rules  described  in  subsection  7.1.  We 
point  out  again  that  we  did  not  know  in  advance  values  refractive  indexes  of  our  wooden 
cubes.  Therefore,  when  applying  stopping  rules,  we  were  unbiased.  Table  3  presents  numbers 
an  and  bn  =  an/an_i  for  the  case  1.1(1)  (see  Table  2  for  labeling  of  our  cases).  Behavior 
of  these  numbers  for  other  cases  was  similar.  Therefore,  Table  4  presents  only  numbers 
for  the  final  iteration.  Figures  6  and  7  display  computed  images.  Our  goal  was  twofold:  (1) 
to  obtain  accurate  locations  of  inclusions  and  (2)  to  accurately  image  refractive  indexes  in 
them.  However,  we  did  not  have  a  goal  to  accurately  image  shapes  of  inclusions. 

Table  3.  Computational  results  for  the  case  1.1(1),  see  Table  2  for  labeling  of  cases  and  (49)  for  the  stopping 
rule 


Iter.,  n 

-hO 

or 

Ctn 

bn 

£f 

nf  =  V 'el 

2 

1.28 

0.027 

0.21 

3 

2.53 

0.209 

7.74 

4 

2.9 

0.160 

0.76 

5 

3.76 

0.266 

1.66 

6 

4.66 

0.580 

2.18 

=  £{r]  =  4.66 

2.16 

7 

5.6 

0.683 

1.18 

8 

8.1 

0.809 

1.18 

7.3  Accuracy  of  the  blind  quantifiable  imaging 

We  have  independently  measured  refractive  indexes  after  the  above  images  were  obtained. 
Those  measurements  were  performed  by  two  well  established  methods:  the  Waveguide 
Method  [19]  and  the  Oscilloscope  Method,  see  [12]  and  http://www.tek.com/learning/oscilloscopes 
for  the  latter  method.  Tables  5  and  6  display  comparisons  of  our  computational  results  with 
measured  ones,  along  with  both  measurement  and  imaging  errors.  One  can  see  the  error  in 
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Table  4.  Computational  results  for  five  cases,  see  (48)  for  the  stopping  rule  and  Table  2  for  labeling  of  cases. 
The  rest  of  iterations  for  all  these  five  cases  was  similar  with  Table  3.  Comparison  of  Table  4  with  (48) 
makes  it  clear  why  either  of  function  £r™')  or  e(”  1'1  was  chosen  as  the  final  imaging  result  e( 


our  computations  did  not  exceed  the  measurement  error  in  all  cases  except  of  the  case  1.1(2) 
in  Table  6,  in  which  our  error  was  1.8%  more  than  the  measurement  error. 

Table  5.  Comparison  of  imaging  results  of  values  of  refractive  indexes  for  six  cases  of  Table  2  with  measure¬ 
ments  by  the  Waveguide  Method 


Case 

Blindly  imaged  n  :=  rif 

Measured  n,  error 

Imaging  error 

1.1(1) 

2.16 

2.07, 11% 

4.3% 

1.1(2) 

2 

2.07, 11% 

3.4% 

1.2 

2.16 

2.07, 11% 

4.3% 

1.3 

2.19 

2.07, 11% 

5.8% 

2.1 

1.73 

1.71,3.5% 

1.2% 

2.2 

1.79 

1.71,3.5% 

4.7% 

Table  6.  Comparison  of  imaging  results  of  values  of  refractive  indexes  for  six  cases  of  Table  2  with  measure¬ 
ments  by  the  Oscilloscope  Method 


Case 

Blindly  imaged  n  :=  rif 

Measured  n,  error 

Imaging  error 

1.1(1) 

2.16 

2.17,6% 

0.5% 

1.1(2) 

2 

2.17,6% 

7.8% 

1.2 

2.16 

2.17,6% 

0.5% 

1.3 

2.19 

2.17,6% 

1% 

2.1 

1.73 

1.78,6% 

2.8% 

2.2 

1.79 

1.78,6% 

0.56% 

8  Discussion 

The  goal  of  this  research  effort  was  to  verify  the  globally  convergent  numerical  method  of 
[2,  3,  4,  5]  experimentally.  We  have  arranged  picosecond  scale  time  resolved  measurements 
to  acquire  the  data  scattered  by  small  dielectric  abnormalities  (lps  =  10_12sec  =  10~3ns). 


a)  maxn^1)  =  1.02 


b)  maxt^2'  =  1.13 


c)  maxrh3)  =  1.59 


d)  maxn^  =  1.70 


e)  ma xn^  =  1.94  f)  maxn^6^  =  max  nf  =  2 

Figures  6-a)-h)  show  the  dynamics  of  the  sequence 
ues  of  refractive  indexes  maxpn^^  =  \J maxp  ei^ 
{x:nSk\x)  =  maxpn^k\x)}.  The  final  image 
plodes”  on  the  second  iteration  after  the  stop,  st 


of  images  for  the  case  number  1.1(1).  Maximal  val- 
are  displayed.  Each  image  represents  the  level  surface 
is  presented  on  f).  h)  shows  that  the  image  “ex- 
se  the  Stopping  rule  (48)  and  Table  3. 


The  time  resolved  data  were  measured  in  the  tomographic  manner  on  a  planar  surface,  which 
was  opposite  to  the  single  source  location  used.  This  constitutes  to  a  minimal  amount  of 
information  gathered,  which  is  important  for  such  applications  as,  e.g.  imaging  of  land  mines, 
imaging  of  baggages  in  airports,  more  generally,  imaging  of  explosives  and  also  imaging  of 
defects  in  non-distractive  testing.  Furthermore,  the  technique  of  [2,  3,  4,  5]  can  be  extended 
to  the  most  practically  interesting  case  of  stand  off  detection  using  the  back  reflected  data 
only  [3]. 

At  each  location  of  the  detector  the  data  were  collected  for  total  of  12,300  picoseconds, 
i.e.  12.3  nanoseconds.  A  radically  new  data  processing  procedure  was  developed.  The  goal 
of  this  procedure  was  to  adapt  the  data  to  the  mathematical  model,  which  is  based  on  the 
hyperbolic  equation  (1).  Results  of  the  data  processing  were  used  as  in  input  for  the  globally 
convergent  numerical  method  of  [2,  3,  4,  5]  for  the  CIP  (1),  (2),  (6).  While  processing 
the  data,  we  have  used  the  fact  that  on  five  out  of  six  sides  of  the  prism  0  the  data  were 
not  sensitive  to  the  presence  of  inclusions  and  thus  we  have  assigned  values  taken  from 
the  reference  medium  for  them.  The  only  sensitive  side  was  the  planar  surface  P  of  our 
measurements.  However,  as  soon  as  the  data  are  pre-processed,  our  algorithm  does  not  use 
any  a  priori  knowledge  neither  about  the  background  values  of  the  unknown  coefficient,  nor 
about  locations  of  inclusions.  The  only  information  it  uses  is  that  the  value  of  the  refractive 
index  outside  of  the  computational  domain  of  interest  is  the  same  as  one  in  the  air  and  that 
the  refractive  index  in  the  computational  domain  is  not  less  than  the  one  in  the  air.  Our 
algorithm  solves  a  fully  nonlinear  problem.  No  linearization  is  used.  Our  computations  were 
quite  fast  and  took  only  a  few  minutes  on  a  regular  PC  to  get  the  above  images.  This  points 


a)  Case  1.1(1),  maxnj  =  2.16  b)  Case  1.1(2),  ma xn/  =  2  c)  Case  1.2,  ma xn/  =  2.16 
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d)  Case  1.3,  ma xnf  =  2.19  e)  Case  2.1,  maxnj  =  1.73  f)  Case  2.2,  maxnj  =  1.79 

Figures  7-a)-f)  display  resulting  images.  It  should  be  kept  in  mind  that  we  did  not  have  a  goal  to  image  shapes 
of  inclusions  accurately.  Rather,  our  goal  was  only  to  image  their  locations  and  maximal  values  of  refractive 
indexes  n/(x)  =  \fel .  On  each  figure  n/[x)  =  maxnj  for  all  points  of  the  image  of  the  corresponding 
cube.  In  addition  to  the  cut-offs  (47),  we  have  made  the  last  post  processing  cut-off  of  the  imaged  function 
on  each  figure  just  to  make  it  look  better.  That  cut-off  was  made  around  the  center  of  the  image.  For  all 
cases  the  dynamics  of  the  change  of  images  of  functions  with  iterations  was  similar  with  one  on  Figures 
6-a)-h). 

towards  a  possibility  to  obtain  real  time  imaging  after  a  certain  programming  effort. 

We  have  formulated  the  global  convergence  theorem  of  [3]  for  our  method  and  have 
discussed  its  features  in  detail.  This  theorem  has  a  problem  with  the  truncation  of  the  pseudo 
frequency  s,  since  we  cannot  prove  convergence  as  s  — >  oo.  Furthermore,  if  we  would  somehow 
prove  this  convergence,  then  we  would  address  a  long  standing  and  well  known  question 
about  uniqueness  theorem  for  the  Coefficient  Inverse  Problem  (1),  (2),  (6).  In  addition,  we 
have  pointed  out  that  a  similar  issue  exists  in  the  classic  Real  Analysis  in  connection  with 
asymptotic  series  expansions  of  some  functions,  although  still  those  expansions  are  known 
to  provide  good  approximations.  We  have  also  pointed  out  to  a  similar  problem  with  the 
truncation  of  the  Fourier  series  in  the  Gelfand-Krein-Levitan  reconstruction  method  for  a 
2-D  hyperbolic  CIP,  although  numerical  results  demonstrate  an  excellent  performance  of  this 
method  [14]. 

Thus,  we  have  concluded  that  the  only  way  to  justify  such  truncations  is  via  numerical 
studies.  This  was  done  successfully  in  [2,  3,  4,  5]  for  computationally  simulated  data  and  also 
for  the  experimental  data  in  the  current  paper.  We  believe  that  the  success  in  these  numerical 
studies  has  the  same  nature  as  routine  truncations  of  high  frequencies  in  engineering.  Indeed, 
it  is  well  known  that  in  engineering  high  frequencies  are  often  truncated  without  any  proofs 
of  convergence  and  things  still  work  well. 

Since  the  question  about  convergence  at  s  — >  oo  is  very  challenging  to  address,  we  have 
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proposed  a  new  mathematical  model  for  our  numerical  method.  We  believe  that  this  model 
intrinsically  has  the  same  nature  as  the  above  mentioned  truncation  of  asymptotic  series. 
By  our  model,  a  certain  small  positive  parameter  £  is  allowed  to  be  infinitely  small  inde¬ 
pendently  on  the  truncation  pseudo  frequency  s,  although  formally  it  does  depend  on  s. 
Furthermore,  it  turns  out  that  we  have  benefited  from  this  model  via  applying  it  in  compu¬ 
tations  (subsection  6.1).  It  is  also  important  to  point  out  that  all  features  of  our  numerical 
implementation  mentioned  in  subsection  6.1  were  pre-arranged  prior  getting  experimental 
data.  Therefore,  we  were  not  biased  when  computing  above  results.  Our  stopping  criterion, 
which  truncates  iterations,  is  in  a  full  agreement  with  one  of  backbone  ideas  of  the  theory 
of  Ill-Posed  Problems,  see  pages  156  and  157  in  [9].  By  this  idea,  the  iteration  number  can 
be  used  as  a  regularization  parameter. 

In  fact,  we  had  semi  blind  testing.  In  other  words,  we  knew  locations  of  inclusions, 
although  this  information  is  not  used  in  our  algorithm.  However,  we  were  unaware  about 
values  of  refractive  indexes  in  them.  So,  with  respect  to  the  values  of  the  refractive  indexes  we 
were  blind.  We  have  verified  a  posteriori  results  of  our  computations  via  direct  measurements 
using  two  well  established  methods  of  Physics.  These  measurements  have  revealed  that  our 
computational  results  for  experimental  data  have  consistently  demonstrated  an  excellent 
accuracy  in  six  (6)  out  of  six  (6)  cases,  i.e.  in  100%  available  cases.  In  all  cases,  our 
computational  error  was  either  less  or  no  more  than  1.8%  that  the  measurement  error.  An 
accurate  image  was  obtained  even  in  the  most  difficult  case  1.3  in  which  the  inclusion  was 
located  in  an  area  where  the  amplitude  of  the  EM  wave  was  much  less  compared  with  the 
one  on  the  center  line. 

Actually,  we  had  quite  a  large  noise  in  our  input  data  for  the  globally  convergent  method. 
This  noise  was  generated  by  the  following  four  factors: 

1.  Our  governing  PDE  (1)  cannot  be  derived  from  the  Maxwell  system. 

2.  Our  theory  does  not  work  for  our  case  when  the  coefficient  er  has  a  discontinuity 
at  the  inclusion/background  interface.  Items  1  and  2  indicate  that  we  have  worked  with  a 
simplified  mathematical  model. 

3.  The  experimental  data  itself  naturally  have  a  large  noise. 

4.  We  had  an  implicit  noise  in  our  follow  up  mathematical  model  occurred  in  the  data 
pre-processing. 

Thus,  we  believe  that  the  excellent  accuracy  of  our  results  of  blind  testint  points  towards 
a  high  degree  of  robustness  of  the  technique  of  [2,  3].  In  summary,  based  on  above  results,  we 
conclude  that  the  globally  convergent  numerical  method  of  [2,  3]  is  now  completely  validated. 

9  Ongoing  Effort 

After  the  above  results  were  obtained,  we  have  undertaken  an  additional  effort  with  the 
goal  to  refine  these  results,  especially  to  refine  shapes  of  imaged  inclusions.  This  is  a  time 
consuming  effort,  as  any  effort  of  a  new  research  software  development.  So,  the  effort  outlined 
in  this  section  is  currently  ongoing. 
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9.1  The  two-stage  numerical  procedure  of  [3,  4,  5] 


In  previous  works,  which  were  submitted  for  publication  in  2009  [3,  4,  5],  a  two-stage  nu¬ 
merical  procedure  was  developed.  On  the  first  stage,  the  globally  convergent  method  of  [2] 
has  provided  a  good  first  approximation  for  the  solution.  And  on  the  second  stage,  this 
approximation  was  taken  as  the  starting  point  for  further  refinement  by  a  locally  convergent 
technique.  The  key  point  here  is  that  there  is  a  rigorous  guarantee  that  this  first  approx¬ 
imation  is  close  to  the  exact  solution.  Another  important  point  is  that  not  every  locally 
convergent  technique  works  here.  Namely,  it  was  established  in  numerical  experiments  of 
[3,  4,  5]  that  the  application  of  the  quasi-Newton  method  on  the  same  mesh  as  one  used  on 
the  globally  convergent  part  does  not  lead  to  a  refinement  of  the  image.  Thus,  in  [3,  4,  5] 
the  Adaptive  Finite  Element  technique  (adaptivity)  was  used.  Namely,  mesh  was  refined 
locally  on  several  iterations.  And  then  the  quasi-Newton  method  was  applied  on  this  refined 
mesh.  Rigorous  criteria  for  mesh  refinements  were  established  in  [4,  5].  Specifically,  mesh 
should  be  refined  at  regions  where  the  values  of  the  modulus  of  the  gradient  of  the  Tikhonov 
regularization  functional  are  close  to  maximal  ones. 

So,  the  procedure  of  [3,  4,  5]  was  recently  applied  to  the  above  experimental  data  [1]. 
While  on  Figures  7  we  have  reconstructed  only  locations  of  inclusions  and  refractive  indexes 
in  them,  now  we  can  reconstruct  shapes  of  inclusions  as  well,  see  Figures  8.  Thus,  we 
now  have  everything  we  want  from  experimental  data:  (a)  locations  of  abnormalities,  (b) 
refractive  indexes  in  them,  and  (c)  their  approximate  shapes.  Of  course,  it  is  hard  to  expect 
that  more  sophisticated  shapes  will  also  be  reconstructed.  But  at  least,  we  anticipate  that 
we  will  be  able  to  reconstruct  convex  hulls  of  abnormalities,  i.e.  we  would  approximate  their 
sizes. 


Figure  8.  A  sample  of  the  image  from  experimental  data  via  applying  the  two-stage  procedure. 

An  important  point  to  make  here  is  that  while  using  the  adaptivity,  we  have  used  only 
the  data  in  time  domain  without  calculating  Laplace  transform.  The  procedure,  similar  with 
one  depicted  on  Figure  6  (subsection  5.3.2)  was  used  in  time  domain,  for  each  time  step. 
The  second  important  point  is  that  while  using  the  adaptivity,  we  need  to  rely  on  a  priori 
known  upper  estimate  for  the  unknown  coefficient  er  (a;)  :  otherwise,  imaged  values  of  this 
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coefficient  would  be  higher  than  they  should  be.  Thus,  we  took  those  upper  estimates  from 
the  globally  convergent  part. 


9.2  The  quasi-reversibility  method  for  global  convergence 

The  idea  is  to  work  with  the  boundary  data  which  are  given  at  a  part  of  the  boundary  only, 
rather  than  on  the  entire  boundary.  The  above  described  version  of  the  globally  convergent 
method  works  for  the  case  when  one  has  forward  scattering  data  (Figure  1)  and  the  Dirichlet 
boundary  condition  given  at  the  entire  boundary.  Although,  for  the  case  of  our  experimental 
data,  we  have  managed  to  assign  the  Dirichlet  boundary  condition  from  the  reference  medium 
at  the  rest  of  the  boundary,  it  is  interesting  to  investigate  the  case  when  only  the  radiation 
boundary  condition  for  the  function  w  in  (8)  is  given  at  the  rest  of  the  boundary.  One  can 
prove  that  this  condition  is  transformed  in  the  Neumann  boundary  condition  for  the  function 
Q- 

On  the  other  hand,  the  above  described  experiment  has  provided  us  only  with  the  Dirich¬ 
let  boundary  data  at  the  bottom  side  P  of  the  domain  O.  Let  P  =  {z  =  —2.4}  be  the  plane 
containing  the  rectangle  P.  Having  the  Dirichlet  boundary  data  at  P  and  extending  them 
as  one  for  the  reference  medium  on  the  rest  of  P,  one  obtains  the  Dirichlet  boundary  data 
Dp  (x,y)  at  the  entire  plane  P.  Then  one  can  solve  the  following  boundary  value  problem 
in  the  half  space  below  P 


Aw  —  s2w  =  0,  (x,  y,z)  G  { (x,  y )  G  M2,  z  <  —2.4}  , 
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w  | p=  Dp  (x,y). 


The  solution  of  this  problem  provides  one  with  the  Neumann  boundary  condition  at  P  for 
the  function  w.  In  particular,  one  obtains  at  the  rectangle  P  C  P 


Wz  | P=  m  (x,y,  s) , 


where  the  function  m  is  known. 

This  means  that  we  have  obtained  both  Dirichlet  and  Neumann  boundary  conditions 
for  the  function  q  at  the  bottom  side  P  of  the  domain  D.  It  is  very  important  that  similar 
considerations  can  be  applied  to  the  most  interesting  case  of  back  reflected  data  [3] . 

On  the  other  hand,  the  quasi-reversibility  method  [15,  17]  is  well  suited  for  solving  PDEs 
with  over-determined  boundary  conditions  on  one  part  of  the  boundary  and  no  boundary 
conditions  on  the  rest  of  the  boundary  [15,  17].  Another  thought  behind  the  application  of 
this  method  is  that  it  over-smooths  solutions.  So,  this  would  be  an  alternative  to  the  above 
mentioned  adaptivity  procedure  in  our  goal  to  improve  shapes  of  imaged  inclusions. 

Let  Ln  be  the  operator  in  the  left  hand  side  of  (17)  and  Fn  be  the  right  hand  side  of  (17), 
in  which  the  nonlinear  term  (Vgn)2  is  properly  replaced  as  described  in  subsection  3.1.  Then 
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in  the  quasi-reversibility  method  one  minimizes  the  following  Tikhonov  functional  [15,  3,  4] 

Jot  (?n)  —  II Ln  (<?n)  —  ^n||x,2(n)  +  OL  ||?n||#2(Q)  , 

subject  to  boundary  conditions,  ffere  a  is  the  Tikhonov  regularization  parameter  [24],  To 
minimize  this  functional,  finite  differences  are  used  and  the  minimization  is  carried  out  over 
values  of  the  function  qn  at  grid  points.  When  working  with  back  reflected  data,  it  was 
observed  in  [3]  that  one  should  use  the  above  described  two-stage  procedure.  On  the  first 
stage  the  globally  convergent  method  provides  good  locations  of  inclusions.  And  on  the 
second  stage  the  gradient  method  provides  accurate  values  of  the  coefficient  er  (x)  within 
these  inclusions.  Figure  9  represents  the  result  obtained  in  [3]  for  the  case  of  back- reflected 
data.  The  plane  wave  falls  from  the  top  and  measurements  of  the  back-reflected  signal  are 
carried  out  on  the  top  boundary.  The  radiation  boundary  condition  is  imposed  at  other 
parts  of  the  boundary,  which  transforms  into  the  Neumann  boundary  condition  for  functions 

qn- 


o 


Figure  9:  The  application  of  the  quasi-reversibility  to  the  case  of  back  reflected  data  collected  at  the  top  side 
of  the  boundary  (computational  simulations,  see  the  text  for  some  details).  We  have  modeled  imaging  of 
antipersonnel  land  mines,  whose  depths  does  not  exceed  10  centimeters.  Fig.  9a  displays  the  correct  image 
of  two  mine-like  targets  with  er  =  6  in  the  left  target  and  er  =  4  in  the  right  one  and  er  =  1  in  the  rest 
of  the  medium.  Fig.  9b  shows  the  image  obtained  after  applying  the  globally  convergent  method  with  the 
quasi-reversibility.  Locations  of  both  mines  are  imaged  accurately.  However,  values  of  er  are  about  60% 
of  the  correct  ones.  Figure  9c  displays  the  final  image  obtained  on  the  second  stage  of  the  procedure.  Both 
locations  of  targets  and  values  of  er  inside  of  them  are  accurately  imaged,  so  as  the  value  er  =  1  outside  of 
these  targets. 


10  Presentations 

Results  of  our  work  on  experimental  data  [2]  were  presented  by  the  PI  at  six  (6)  professional 
meetings  listed  below.  While  titles  of  presentations  were  varied  a  somewhat,  the  commonly 
used  title  can  be  labeled  as  “Picosecond  scale  experimental  verification  of  a  globally  conver¬ 
gent  algorithm  for  a  coefficient  inverse  problem” . 
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Specifically,  these  results  were  presented  at: 

1.  Army  Research  Laboratory 

2.  Seminars  of  the  Department  of  Mathematics  of  the  University  of  North  Carolina  at 
Charlotte 

3.  SPIE  meeting  in  San  Francisco 

4.  Henry  Poincare  Institute,  Paris,  France 

5.  Workshop  on  Inverse  Problems  at  Chinese  University  of  Hong  Kong  (plenary  speaker) 

6.  International  Conference  on  Inverse  Problems,  Antalia,  Turkey 

11  Honor 

The  publication  [2]  was  selected  as  a  featured  paper  of  Inverse  Problems,  because  of  a  high 
rating  by  referees. 
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